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EXECUTIVE  SUMMARY 


The  pressure  modeling  technique  is  used  to  study  upward  fire  spread  on  fuel 
walls  composed  of  char-forming  or  laminated  materials.  Time-resolved 
measurements  are  obtained  at  one-atmosphere  (full-scale)  and  at  elevated 
pressure  (model  scales)  to  characterize  fire  growth  in  terms  of  rate  of  total 
mass  loss,  flame  height  and  maximum  lateral  flame  extent  during  the  spread 
process.  The  char-forming  materials  (pine-wood,  particle-board  and  a  rigid, 
polyurethane  foam)  are  tested  in  a  90°  wall  corner  configuration  while  the 
laminated  materials  (polymethylmethacrylate  (PMMA)  or  ceramic  backings)  are 
tested  in  a  wall  configuration.  Thermally- thick  PMMA  is  tested  in  both 
configurations  for  purposes  of  comparison. 

The  following  results  are  obtained  from  the  modeling  study: 

1.  Pressure  modeling  is  sufficiently  accurate  for  the  prediction  of  fire 
growth  from  a  point  ignition  on  a  uniform  PMMA  wall  when  both  upward  and 
lateral  flame  spread  processes  occur. 

2.  The  behavior  of  the  flame  spread  process  at  elevated  air  pressures,  for 
walls  composed  of  a  face  layer  of  PMMA  with  a  thick  nonburning  backing 
layer,  is  not  completely  consistent  with  a  simplified  analysis  of  thermal 
conduction  processes. 

3.  Pressure  modeling  of  fire  growth  in  a  wall  corner  configuration  is 
quite  accurate  for  cellulosic,  char-forming  materials  and  for  PMMA.  The 
cellulosics  at  one  atmosphere  exhibit  a  flame  extinction  phenomenon  due  to 
char  buildup  after  significant  lateral  flame  spread  that  is  not  observed 
at  elevated  pressures. 

4.  Thermally  thick,  rigid,  high  density  polyurethane  foam  in  a  corner 
configuration  will  not  support  significant  flame  spread  at  one-atmosphere 
but  will  at  elevated  pressures  with  a  properly  scaled,  small  PMMA  ignition 
source.  This  behavior  is  perhaps  due  to  excessive  radiant  heat  loss  from 
the  char  and  the  intumescent  character  of  the  char  at  one-atmosphere.  Gas 
phase  chemical  kinetics,  which  may  be  the  most  important  factor  in  the 
initiation  of  flame  spread  on  the  full-scale  foam,  is  clearly  not 


modeled . 


5.  A  simplified  analysis  of  thermal  conduction  in  a  laminated  material  is 
used  to  show  how  flame  spread  rates  are  affected  by  the  thermal  properties 
of  the  face  layer  and  backing  at  both  one- at mo sphere  and  at  elevated 
pressures. 

6.  A  numerical  technique  is  used  to  predict  one-dimensional,  transient 
fuel  mass  flux,  fuel  surface  temperature  and  char  thickness  during 
material  exposure  to  a  prescribed  radiant  (and  convective)  heat  flux. 
Calculated  results  show  that  reasonable  pressure  modeling  of  flame  spread 
rates  should  be  expected  for  the  cellulosic  fuels  due  to  increases  in 
surface  temperature  and  char  production  for  conditions  simulating  elevated 
air  pressure. 

These  results  lead  to  two  main  conclusions: 

1.  Pressure  modeling  of  three-dimensional  fire  spread  on  uniform  wal 
and  wall  corners  has  now  been  validated  for  PMMA  and  for  wood  fuels, 
has  not  yet  been  established  that  the  modeling  technique  is  valid  for 
predicting  fire  growth  on  other  charring  fuels  in  corner  configurations. 
Results  from  this  and  previous  studies  have  shown  that  in  general,  the 
complex  process  by  which  self-sustained  fire  spread  is  initiated  is  not 
pressure  modeled.  Such  fire  spread  initiation  occurs  readily  at  elevated 
pressures  because  surface  radiant  heat  loss  and  the  action  of  gas  phase 
chemical  retardents  cannot  be  modeled.  With  the  wood  and  PMMA  fuels,  fire 
spread  rates  on  wall-ceiling  corners  should  also  be  predictable  by 
pressure  modeling,  based  on  previous  work^^  with  ceiling  channel 
configurations . 

2.  It  appears  that  much  more  experimental  information  is  needed  before 
pressure  modeling  can  be  used  to  predict  fire  growth  on  practical 
composite  materials.  At  present  the  thickness  of  all  components 
(including  adhesives)  within  the  thermal  wave  developed  during  fire  spread 
must  be  modified  according  to  the  pressure  modeling  scheme.  However, 
radiant  exposure  in  real  enclosure  fires  may  well  be  sufficiently  high  (2- 
U  W/cm  )  to  confine  thermal  wave  penetration  to  a  surface  layer  of  fuel 
during  fire  spread.  Pressure  modeling  of  such  a  fire  spread  process  would 
then  be  accurate  without  any  modification  of  the  laminated  material. 
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INTRODUCTION 


1.1  PURPOSE 


The  purpose  of  this  project  is:  1)  to  study,  by  experiment,  the  behavior  of 
upward  fire  spread  and  growth  at  both  atmospheric  and  elevated  air  pressures; 
and  2)  to  perform  an  analysis  of  transient,  one-dimensional  heat  conduction 
and  pyrolysis,  in  order  to  determine  the  applicability  of  the  pressure  model¬ 
ing  technique  to  char-forming  and  laminated  materials. 

1.2  BACKGROUND 


Materials  within  an  aircraft  cabin  can  be  exposed  to  flame  heat  transfer  from 
external  fuel-spill  fires  after  an  aircraft  accident.  It  is  desirable  to  have 
cabin  materials  which  will  limit  any  compartment  fire  growth  initiated  by  the 
external  fire.  Methods  for  characterizing  the  fire  growth  potential  of  a 
variety  of  aircraft  cabin  materials  in  various  configurations  are  therefore 
needed.  One  method  for  studying  the  effect  of  fuel  configuration  and  geometry 
on  transient  fire  growth  is  the  pressure  modeling  technique,  whereby  experi¬ 
ments  with  small-scale  fuels  are  conducted  at  elevated  air  pressures  to 
simulate  full-scale,  controlling  fluid  mechanic  and  thermal  parameters.  It  is 
important  to  determine  how  effects  due  to  changes  in  fuel  composition  will 
alter  such  a  modeling  process. 

Fire  growth  within  an  aircraft  cabin  can  occur  by  several  different  modes. 

One  such  mode  is  by  upward  and  lateral  fire  spread  on  vertical  cabin  sur¬ 
faces.  In  a  previous  study^ ^ ,  the  feasibility  of  modeling  fire  spread  on 
vertical  walls  through  the  use  of  small-scale  experimi  nts  at  elevated  air 
pressure  (pressure  modeling  technique)  was  proven  for  a  homogeneous  fuel  with 
simple  vaporization  characteristics,  polymethylmethacrylate  (PMMA).  The 
present  study  deals  with  fuels  which  undergo  pyrolysis  reactions  leading  to 
char  formation  and  fuels  that  are  composites  of  two  different  materials  lamin¬ 
ated  together.  While  such  materials  are  more  like  real  fuels  than  homogeneous 
PMMA,  an  effort  has  been  made  to  use  the  simplest  possible  charring  and  lamin¬ 
ated  materials  so  that  the  pressure  modeling  behavior  can  be  understood. 
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Relatively  simple  fuel  configurations  are  also  used  for  this  study  in  order  to 
facilitate  subsequent  analysis.  Laminated  materials  are  burned  in  an  unbound¬ 
ed  wall  configuration  ignited  at  a  single  point  to  take  advantage  of  the  flame 
spread  characteristics  of  the  surface  layer  of  PMMA  fuel.  Char-forming 
materials  will  ordinarily  not  sustain  extensive  flame  spread  unless  exposed  to 
a  minimum  heat  flux  level.  In  the  present  study,  the  exposure  is  provided,  in 
part,  by  a  small,  PMMA  initiator  built  into  the  charring  fuel.  The  remainder 
of  the  exposure  flux  results  from  the  use  of  a  90°  corner  configuration,  which 
allows  some  thermal  radiation  lost  from  one  wall  of  the  corner  to  be  captured 
by  the  opposite  wall. 


PRESSURE  MODELING  EXPERIMENTS 

2.1  EXPERIMENTAL  ARRANGEMENT  SKETCHED  IN  FIGURES  1  AND  3 

2.1.1  CONFIGURATION  OF  LAMINATED  FUELS.  An  unbounded  vertical  wall  ignited 
at  a  point  near  the  wall  base  is  used  to  test  pressure  modeling  concepts  for 
laminated  fuels  and  for  uniform  PMMA.  Tables  1  and  2  list  the  fuel  composi¬ 
tions  and  dimensions,  respectively,  used  in  the  assembly  of  the  full-scale  and 
model  wall  configurations.  Those  dimensions  in  Table  2  denoted  as  "scaled" 
(fuel  height  and  width,  etc.)  are  reduced  from  the  full-scale  values  listed 
when  used  in  the  model  configurations.  As  explained  in  References  1  and  2, 
the  pressure  modeling  scheme  requires  that  such  a  reduction  of  length  scales 
be  made  proportional  to  the  minus  2/3  power  of  absolute  air  pressure.  At  air 
pressures  of  11.2,  20.5  and  31.6  atmospheres,  this  tranlates  into  respective 
length  scale  reductions  by  factors  of  1/5,  1/7.5  and  1/10.  Fuel  thickness, 
except  for  that  of  the  uniform  PMMA  wall,  is  not  scaled  in  this  manner  due  to 
practical  limitations.  Instead,  a  thermally  thick  backing  layer  is  used  while 
the  PMMA  face  layer  is  maintained  at  a  constant  thickness  for  all  pressures. 

A  fixed  fuel  thickness  would  be  used  in  practice  if  complex  composite  fuels 
were  to  be  subjected  to  any  flammability  test.  For  the  case  of  uniform  PMMA, 
cast  slabs  1/5  and  1/10  the  full-scale,  31.75  nm  thickness  happen  to  be 
commercially  available.  Another  exception  to  the  modeling  scheme  is  the  case 
of  the  PMMA- PMMA  laminate  wall.  The  model  fuels  correspond  to  the  PMMA-Inert 
prototype  dimensions,  instead  of  those  for  the  PMMA-PMMA  laminate. 
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TABLE  I 

DESCRIPTION  OF  WALL  MATERIALS  AT  AMBIENT  TEMPERATURE 


Wall  Material 

Composition 

Density 

[kg/m3! 

Thermal  Diffusivity 
[m2/s J 

PCX 

[J2/mVs] 

Inert 

"Cotronlcs"  type 

360  Ceramic  Board 

256 

2.425  x  10-7 

1.742  x  10' 

PMMA 

(Melting) 

Polymethylmethacrylate 

("Plexiglas,  cast, 

type  G) 

1200 

7.95  x  10-8 

5.54  x  105 

Pine-Wood 

(Char-Forming) 

Cellulose  "1  x  10" 

pine  board 

434 

1.92  x  10*7 

5.21  x  104 

Particle-Board 

(Char-Forming) 

Cellulose  particleboard 

( "Versaboard"  Douglas  Fir) 

694 

9.07  x  10'8 

1.75  x  105 

CM- 3  7 

(Char-Forming) 

Rigid  Polyurethane  Foam 

(from  NBS  office  of 

Standard  Reference 

Materials)  Polymeric 

Isocyanate 

320 

1.76  x  10-7 

1.98  x  104 
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TABLE  II 


Fuel 

Typg _ 

PMMA  Wall 

PMMA-PMMA 
laminate  wall 

PMMA- Inert 
laminate  wall 

PMMA  wall 
corner 

Particle¬ 
board  wall 
corner 

Pine-wood 
wall  corner 

Inert  wall 
corner 

GM-37 

wall  corner 


DIMENSIONS  OF  FUEL  CONFIGURATIONS 
SKETCHED  IN  FIGURES  1  AND  5 

Scaled  Scaled  Width  of  Wall 

Material  Thickness  Height  (each  leg  of  corner) 

_ [mail _ (mj _ jmj _ 

31.75,  scaled  0.8985  0.2286 

3.18  -  front  face  0.8985  0.2286 

31.75  -  backing  at  1  atm 
12.7  -  backing  at  >  1  atm 

3.18  -  front  face  1.22  0.305 

25. A  -  inert  backing 


31.75  at  1  atm  0.8985  0.137 

12.7  at  >  1  atm 


19.05  1.22  0.2  and  0.15 


19.05  1.21  0.22 


12.7  0.61  0.2 


51  at  1  atm  1.22  0.2 

25  at  >  1  atm 


PMMA  Initiator 
Scaled  Height 
_ [mm  j 


102 


51 


102 


102 


4 


Both  laminated  fuel  walls  are  assembled  with  a  ceramic  adhesive  ("Cotronics" 
type  9011  to  bond  the  face  layer  to  the  backing  layer.  This  bonding  is  aided 
by  mechanical  fasteners  (machine  screws  at  1  atm  and  small  diameter  bolts  at 

11.2  atm)  or  clamps.  For  the  PMMA-PMMA  laminate,  the  ceramic  adhesive  also 
acts  as  an  inerting  agent,  preventing  the  involvement  of  the  PMMA  backing  in 
fuel  vaporization  while  allowing  heat  transmission  to  the  backing. 

Two  difficulties  were  encountered  during  use  of  the  ceramic  paste  in  the  full- 
scale  laminate;  a  nonuniform  paste  thickness  and  an  overly  long  drying  time, 
both  due  to  the  large  area  of  PMMA  involved.  The  thickness  of  the  ceramic 
paste  varied  from  0.8  mm  up  to  3  mm,  with  the  average  thickness  being  1.6  mm 
at  full-scale  but  perhaps  an  order  of  magnitude  less  in  the  model  fuels. 
Complete  drying  of  the  ceramic  paste  was  prevented  by  the  impermeable  PMMA 
surface  at  full-scale  but  reasonable  drying  times  of  several  hours  were 
possible  with  the  model  fuels. 

2.1.2  CONFIGURATION  OF  CHAR-FORMING  FUELS.  An  unbounded,  90°  degree  corner 
arrangement  is  used  to  test  pressure  modeling  predictions  for  char-forming 
fuels  and  for  uniform  PMMA.  Fire  spread  is  initiated  for  the  char-forming 
fuels  by  a  small  PMMA  corner  insert  which  is  flush  with  the  front  surface  and 
bottom  edges  of  the  larger,  charring  material.  Dimensions  and  compositions  of 
all  these  components  are  listed  in  Tables  1  and  2.  The  width  of  each  leg  of 
the  PMMA  corner  initiator  is  one-half  the  scaled  height  listed  in  Table  2 
while  the  initiator  thickness  is  a  constant  12.7  mm  for  both  the  full-scale 
and  model  tests.  This  constant  PMMA  thickness,  which  results  in  the  initiator 
approaching  a  thermally  thin  condition  in  the  latter  stages  of  the  full-scale 
experiments,  is  probably  not  important  in  the  modeling  of  fire  spread.  To 
characterize  the  flame  height  from  the  PMMA  initiator  alone,  an  inert  wall 
corner  was  fabricated  with  the  initiator  insert. 

For  the  laminated  fuels,  the  thickness  dimensions  listed  in  Table  2  are  not 
scaled  down  with  increased  pressure  since  a  thermally  thick  condition  is 
maintained.  The  larger-than-required  thickness  of  the  models  enables  the 
cellulosic  corners  to  be  fabricated  conveniently  with  finishing  nails.  A 
solvent-type  cement  is  used  to  fabricate  both  the  PMMA  fuel  corners  and  PMMA 


initiators.  The  ceramic  adhesive  noted  above  is  used  to  bond  the  two  legs  of 
the  inert  corner.  This  same  adhesive  is  applied  to  the  outer,  top  and  bottom 
edges  of  all  the  fuel  samples,  (wall  as  well  as  corners)  thereby  confining 
flame  spread  to  the  front  face. 

Because  celluloslc  fuels  absorb  moisture,  both  the  pine-wood  and  particle¬ 
board  configurations  are  dried  thoroughly  before  each  experiment  in  an  oven 
set  for  90-100°C.  Values  of  density  for  the  char-forming  fuels  listed  in 
Table  1  are  measured  from  oven-dried  samples.  The  model  cellulosic  fuels, 
after  being  cut  from  the  same  board  as  the  full-scale  sample,  are  dried  and 
then  stored  (no  more  than  2  days)  in  plastic  bags  until  being  placed  in  the 
pressure  vessel.  To  insure  a  dry  atmosphere  in  the  vessel,  room  air  is  purged 
from  the  chamber  before  pressurization  with  air  having  a  dew  point  below  194  K. 
(-79°C) . 

2.1.3  EXPERIMENTAL  MEASUREMENTS  AND  PROCEDURES.  The  flame  spread  experiments 
commence  with  the  ignition  of  a  23.4  mm  long,  6.35  mm  diameter  PMMA  rod  at  one 
atmosphere  or  a  small  wooden  toothpick  at  elevated  pressures.  For  the  wall 
configurations,  this  ignition  source  protrudes  from  the  vertical  centerline  a 
scaled  distance  of  25.4  mm  from  the  wall  base,  normal  to  the  fuel  surface. 

The  same  ignition  source  is  used  for  the  wall-corner  configurations.  Instead 
of  being  in  contact  with  the  fuel,  however,  the  rod  or  toothpick  is  inserted 
horizontally  into  the  apex  of  an  inert,  corner-shaped  slab  of  "Cotronics" 
board  upon  which  the  fuel  corner  rests.  The  ignition  source  is  then  directly 
below  the  bottom  edge  of  the  apex  of  the  PMMA  initiator,  separated  from  the 
PMMA  by  12.7  nxn  at  one-atmosphere  and  by  3  mm  at  elevated  pressures.  At  all 
scales,  the  energy  provided  by  the  burning  rod  or  toothpick  is  probably  close 
to  the  minimum  amount  needed  for  initiation  of  upward  fire  spread. 

Once  flame  spread  up  the  wall  or  corner  begins,  flame  position  and  total  fuel 
mass  are  measured  as  functions  of  time.  The  mass-loss  measurement,  obtained 
from  a  strain  gauge  type  of  load  transducer,  contains  random  fluctuations  of 
up  to  _+  0.5  and  ±0.1  grams  at  one  atmosphere  and  elevated  pressures,  respect¬ 
ively.  These  fluctuations  are  0.2%  and  2X  of  the  respective  total  average 
mass-loss  during  the  experiment.  The  relatively  high  noise  level  of  the  mass 
loss  signal  at  elevated  pressures  is  due  to  the  fast  scan  rates  required 
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during  a  maximum  of  30  seconds  of  digital  data  acquisition,  with  a  consequent 
decrease  in  filtering  efficiency.  Most  of  this  noise  is  eliminated  during 
data  reduction  by  obtaining  mass  loss  rates  from  a  least  squares  regression 
fit  to  the  raw  measurements. 

Flame  position  on  or  above  the  fuel  is  obtained  by  motorized,  33  mm  still 

photography  at  rates  up  to  4  frames  per  second.  Projection  of  the  35  mm 

transparencies  onto  a  ground  glass  screen  alows  flame  height  and  width  to  be 
measured  conveniently.  A  vertical  length  scale,  graduated  in  alternately 
colored  1  centimeter  bands,  appears  next  to  the  fuel  in  the  photographs  to 

permit  resolution  of  flame  position  to  within  0.32  of  total  fuel  height. 

Flame  spread  time  is  obtained  from  a  photographed  digital  clock.  At  one-atmo¬ 
sphere,  the  time  measurement  can  be  resolved  to  within  1  second  about  0.22  of 
the  total  fire  spread  duration.  Time  is  resolvable  to  0.1  second  at  elevated 
pressures,  which  represents  about  0.52  of  total  spread  time  in  most  cases. 

2.2  FIRE  GROWTH  BEHAVIOR 

2.2.1  IGNITION  SOURCE.  The  flame  height  above  the  PMMA  rod  ignition  source 
at  one  atmosphere  is  about  0.1  m.  For  the  corner  configuration,  this  means 
that  the  ignition  flame  nearly  reaches  the  top  edge  of  the  PMMA  corner  initi¬ 
ator.  Flame  height  from  the  burning  toothpick  at  11.2  atmospheres  is  0.0175 
to  0.02  m  high,  or  about  the  required  factor  of  five  smaller  than  the  one- 
atmosphere  flame.  At  higher  pressures,  the  toothpick  flame  height  is  roughly 
the  same  20  mn  value  and  thus  is  somewhat  larger  than  required  for  modeling. 

2.2.2  WALL  FIRES.  The  flame  spreads  upward  on  the  wall  from  the  PMMA  igni¬ 
tion  rod  or  wooden  toothpick  at  a  much  higher  rate  than  that  of  lateral  or 
downward  flame  spread.  A  tear-shaped  flame  on  the  surface  of  the  PMMA  fuel 
results,  with  the  maximum  flame  width  occurring  near  the  ignition  level  during 
most  of  the  spread  process  (as  the  total  spread  width  approaches  .14  -  .16  m, 
lateral  spread  rates  near  the  top  of  the  fuel  become  more  important).  The 
uniform  PMMA  wall  is  transparent,  allowing  both  the  flame  and  pyrolysis  fronts 
to  be  seen  by  photography  from  the  (nonburning)  back  side. 
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For  the  PMMA-PMMA  laminate,  the  ceramic  paste  between  the  thin  face  layer  and 
backing  prevents  the  thick  backing  layer  from  burning.  At  one-atmosphere,  a 
region  without  flame  is  visible  at  the  base  of  the  wall  where  the  face  layer 
has  been  consumed.  This  burnout  region  grows  from  .03  m  to  0.3  m  above  the 
wall  base  as  flame  spread  proceeds  from  half-way  up  the  wall  to  the  top  of  the 
fuel.  For  the  PMMA-Inert  laminate,  upward  flame  spread  on  the  face  layer  is 
much  more  rapid  due  to  the  low  density  backing.  The  face  layer  is  consumed 
near  the  base  of  this  laminate  only  after  flame  spread  to  the  top  of  the  wall 
has  occurred.  However,  about  600  seconds  after  ignition,  burnout  of  the  face 
layer  is  observed  to  a  height  of  0.7  m  above  the  wall  base. 

2.2.3  WALL-CORNER  FIRES.  To  characterize  the  0.1  m  high  PMMA  initiator  of 
fire  spread  in  the  corner  configuration,  flame  height  is  measured  with  inert 
walls  supporting  the  PMMA  insert.  At  one-atmosphere,  peak  flame  height  of 
about  0.41  m  occurs  roughly  540  seconds  after  ignition.  Similar  values  of 
flame  height,  about  0.45  m  if  scaled  to  one-atmosphere  conditions,  are  measur¬ 
ed  during  model  tests  at  11.2  and  20.5  atmospheres  for  corresponding  scaled 
times  after  ignition.  Generally,  flame  spread  on  fuel  walls  of  the  corner 
configuration  occurs  within  200-300  seconds,  when  flames  are  just  slightly 
above  the  0.1  m  high  PMMA  initiator. 

Experiments  with  cellulosic  char-forming  fuels  in  the  corner  configuration 
demonstrate  self-sustained  upward  flame  spread  from  a  PMMA  initiator  if  the 
fuels  have  low  moisture  content.  A  51  mm  high  PMMA  initiator  is  found  to  be 
adequate  for  pine-wood  ignition  but  not  for  the  particle-board  fuel,  which 
requires  the  102  ran  high  PMMA.  In  both  cases,  flame  spreads  to  the  top  of  the 
corner  while  generally  confined  to  the  apex  region.  Lateral  spread  then 
follows  most  rapidly  near  the  top  of  the  fuel  as  a  "V"  shaped  pattern  is 
formed.  In  contrast  to  PMMA  fuel,  flame  spread  does  not  continue  laterally  to 
the  outer  edge  of  the  walls  at  one-atmosphere.  The  one-atmosphere  process  of 
lateral  spread  halts  rather  suddenly  with  pine-wood  and  particle-board  when 
the  char  build-up  in  the  apex  region  extinguishes  flaming  combustion  there. 

In  two  tests  with  particle-board  at  one-atmosphere,  complete  extinction  occurs 
reproducibly  after  about  720  seconds  and  a  maximum  lateral  propagation  on  each 
leg  of  0.127  m  from  the  corner  (0.18  m  total  flame  width  projected  onto  a 
plane  normal  to  the  apex  angle  bisector).  Complete  extinction  occurs  with 


full-scale  pine-wood  after  about  360  seconds  and  maximum  lateral  propagation 
of  0.165  m  from  the  corner  (0.23  m  total  projected  flame  width).  This  extinc¬ 
tion  phenomenon  is  not  observed  with  cellulosic  fuels  at  elevated  pressures 
until  well  after  flame  spread  across  the  entire  fuel  surface  is  complete. 

( 3 ) 

Experiments  with  GM-37,  a  rigid  polyurethane  foam  ,  as  the  wall  corner 
material  show  that  self-sustained  upward  flame  spread  is  not  possible  at  one- 
atmosphere  during  exposure  to  0.05,  0.1  or  0.2  m  high  PMMA  initiators.  With 
the  0.1  m  high  initiator  discussed  before,  flame  spreads  upward  to  a  maximum 
height  of  only  0.61  m.  Further  upward  spread  is  retarded  by  bubbling  of  the 
GM-37  surface  (inturaescent  effect)  and  dripping  of  urethane  fuel  down  over  the 
PMMA  initiator.  By  600  seconds  after  ignition,  this  dripping  partially 
extinguishes  the  PMMA  fire.  The  result  for  GM-37  is  very  different  at  elevat¬ 
ed  pressures  since  rapid  upw. rd  and  lateral  flame  spread  occurs. 

2.3  EXPERIMENTAL  RESULTS 

Time  resolved  measurements  of  total  flame  height,  upward  flame  spread  rate, 
maximum  flame  width  and  fuel  mass  loss  rate  are  shown  in  Figures  1-25.  For 
the  transparent  PMMA  fuel  wall,  measurements  of  pyrolysis  zone  height,  upward 
pyrolysis  spread  rate  and  maximum  pyrolysis  width  are  shown  in  addition  to  the 
flame  height  time-history.  Data  from  full-scale  and  model  experiments  has 
been  correlated  in  Figures  1-25  by  utilizing  the  pressure  modeling  scheme. 
Since  the  time  scale  of  model  experiments  is  reduced  as  the  -4/3  power  of 
ambient  air  pressure,  all  time  measurements  are  multiplied  by  the  4/3  power  of 
the  ratio  of  actual  air  pressure,  p,  to  standard  atmospheric  pressure,  pQ. 
Beyond  this  pressure  correction,  the  time  origin  for  many  elevated  pressure 
experiments  has  been  shifted  to  yield  the  best  correlation  of  data  during  the 
initial  stages  of  upward  fire  spread.  The  time  origin  for  the  full-scale 
experiments  corresponds  roughly  to  ignition  of  the  PMMA  fuel  or  corner  initia¬ 
tor. 

2.3.1  FLAME  HEIGHT.  It  can  easily  be  shown  that  successful  modeling  of 
transient  heat  release  rates,  Q  ,  by  elevated  pressure  experiments  also 
implies  successful  modeling  of  flame  heights.  Consider  the  burning  of  fuel 
vapor  in  a  free  plume  above  a  burning  solid.  Empirical  correlations  from 


r 


Steward^  and  from  You  and  Faeth^  give  the  flame  height,  Xj ,  above  such  a 
fuel  source  as  follows: 


x 


f 


<^>2'5 


where  is  ambient  air  density  and  is  a  constant  for  a  given  fuel  heat  of 

combustion,  stoichiometry  and  ambient  temperature.  In  the  pressure  modeling 

scheme,  heat  flux,  Q"  ,  increases  as  p*'  but  since  fuel  area  is  decreasing  as 

p-^^,  total  heat  release  rate  decreases  as  p_2//^.  Ambient  air  density  at 

constant  ambient  temperature  is  simply  proportional  to  ambient  pressure,  p. 

-or » 

As  a  result,  the  flame  height  correlation  yields  a  decrease  in  x^  as  p  , 
which  is  just  the  dependence  required  by  the  pressure  modeling  technique. 


If  combustion  occurs  in  a  wall  plume  (or  wall  fire,  as  in  the  present  experi¬ 
ments)  instead  of  a  free  plume,  the  following  expression  can  be  derived  from 
the  analysis  of  Steward^^  by  assuming  one-half  the  entrainment  of  the  free 
plume: 


x 


f 


c  (^)2/3 

2  V 


where  C2  is  a  constant  for  a  given  fuel  and  O’  is  the  total  heat  release  rate 

per  unit  width  of  the  wall  configuration.  Since  heat  flux,  Q"  ,  should 

2/3  -2/3 

Increase  a  p  while  pyrolysis  height  Xp,  decreases  as  p  , 

x 

O’  ”  fp  Q"  dx  should  be  independent  of  pressure.  The  wall  fire  flame  height 

-2/3 

therefore  decreases  as  p  ,  in  accord  with  the  modeling  scheme. 


-2/3 

Because  all  characteristic  lengths  are  reduced  as  (p/pQ)  ,  the  correction 
of  flame  height,  x^,  (and  pyrolysis  height,  Xp)  for  pressure  by  the  factor 
shown  in  Figures  1-7A  should  result  in  a  correlation  of  data  for  each  materi¬ 
al.  To  the  extent  that  fuel  thickness  is  Important  during  the  flame  spread 
process,  a  high  degree  of  correlation,  and  hence  modeling  success,  would  not 
be  expected  in  those  cases  where  the  sample  thickness  is  the  same  at  both 
atmospheric  and  ambient  pressures. 
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Modeling  of  flame  height  and  pyrolysis  height  on  the  PMMA  wall  (Figures  1  and 
2)  is  quite  good.  It  appears  that  the  maximum  scaled  flame  height  for  the 
models  is  somewhat  greater  than  that  for  the  full-scale  fuel  even  though  the 
fuel  is  thermally  thick  and  fuel  thickness  is  scaled  properly.  Flame  heights 
greater  than  0.9  m,  it  should  be  noted,  are  above  the  top  edge  of  the  PMMA 
wall. 

Figure  3  shows  the  expected  similarity  between  flame  spread  on  the  laminate 
and  that  on  uniform  PMMA  walls  at  one  atmosphere.  This  flame  spread  similar¬ 
ity  is  probably  due  to  the  small  disturbance  of  PMMA  thermal  characteristics 
by  the  ceramic  cement  and  the  small  amount  of  PMMA  surface  layer  consumed 
during  most  of  the  one-atmosphere  flame  spread  process.  At  extinguishment  of 
the  full-scale  PMMA-PMMA  laminate  fire,  the  3.18  mm  face  layer  is  indeed 
consumed  in  a  triangular  region  only  0.46  m  high  and  0.11  m  wide  at  the  base 
of  the  wall.  The  amount  of  face  layer  consumed  was  apparently  not  enough  to 
reduce  the  wall  burning  rate  and  hence  the  flame  height. 

It  is  not  clear  why  the  scaled  flame  heights  for  the  model  PMMA-PMMA  laminates 
are  much  greater  than  the  values  at  one-atraosphere .  In  fact,  these  model 
flame  heights  are  very  similar  to  those  for  the  PMMA-Inert  laminate  in 
Figure  4.  The  inert  backing  in  the  latter  case  reduces  the  heat  loss  suffi¬ 
ciently  to  cause  an  observable  increase  in  flame  spread  rate  compared  to 
uniform  PMMA,  data  for  which  is  also  shown  in  Figure  4.  Another  unexpected 
result  is  the  good  agreement  between  model  and  full-scale  flame  heights  for 
the  PMMA-Inert  laminate.  At  elevated  pressures,  the  PMMA  face  layer  should 
begin  to  simulate  a  thermally  thick  PMMA  slab  rather  than  the  thermally  thin 
condition  that  roughly  exists  at  one-atmosphere.  It  can  be  seen  in  Figure  4 
that  flame  heights  tend  progressively  toward  those  for  the  uniform  PMMA  wall 
as  ambient  pressure  increases  but  this  shift  is  much  less  than  expected. 

Pressure  modeling  accuracy  for  flame  height  on  the  PMMA  wall-corner  is  excep¬ 
tionally  good,  as  shown  by  the  data  in  Figure  5.  This  excellent  flame  height 
correlation  even  extends  to  the  region  above  the  fuel  surface  at  x^  «  0.9  m. 
For  the  particle-board  corner  (see  Figure  6),  modeling  of  flame  height  is  also 
fairly  good  at  11.2  and  20.3  atmospheres.  Flame  heights  are  clearly  less  than 
expected  (below  the  correlation)  at  the  31.6  atmosphere  pressure.  Note  that 
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FIGURE  4  CORRELATION  OF  FLAME  HEIGHT  ON  PMMA-INERT  LAMINATE  WALLS 
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FLAME  HEIGHT,  x((p/p) 


FIGURE  6  CORRELATION  OF  FLAME  HEIGHT  ON  PARTICLE -BOARD  WALL  CORNERS 
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Che  1  atmosphere  extinction  phenomenon  discussed  earlier  leads  to  a  sharp 
decrease  In  full-scale  flame  height  at  about  700  seconds.  A  similar  flame 
height  behavior  Is  seen  In  Figure  7,  where  results  for  pine-wood  are  correlat¬ 
ed.  In  this  case,  the  modeling  scheme  is  quite  accurate  at  all  three  elevated 
pressures.  The  extinction  of  lateral  fire  spread  at  one  atmoshpere  leads  to  a 
peak,  full-scale  flame  height  somewhat  below  that  predicted  by  the  model 
tests . 

Flame  height  measurements  for  the  GM-37  rigid  foam  material  are  correlated  in 
Figure  7A.  The  lack  of  sustained  upward  flame  spread  at  one  atmosphere  is 
evident  from  the  nearly  constant,  0.6  m  flame  height  on  the  1.22  m  high  fuel 
from  350  to  550  seconds  after  ignition.  Flame  heights  associated  with  the 
PMMA  initiator  alone  set  into  an  inert  corner  are  also  shown  in  Figure  7a. 
There  appears  to  be  only  a  small  contribution  from  the  GM-37  fuel  at  one 
atmosphere,  leading  to  a  peak  flame  height  about  0.2  m  greater  than  that  due 
to  the  PMMA  alone.  On  the  other  hand,  upward  flame  spread  is  seen  to  proceed 
rapidly  up  the  total  fuel  height  of  the  models  at  11.2  and  20.5  atmospheres. 
Flame  heights  for  tne  two  model  scales  are  well-correlated  by  the  variables  in 
Figure  7a.  It  is  possible  that  the  full-scale  flame  heights  would  also  follow 
the  model  correlation  if  a  sufficiently  large  initiator  of  flame  spread  were 
used. 


2.3.2  UPWARD  FLAME  SPREAD  RATE.  For  each  of  the  fuel  configurations,  the 

transient  flame  height  data  in  Figures  1-7  are  fit  with  a  cubic  polynomial 

least  squares  regression.  This  polynomial  fit  is  then  differentiated  to  give 

vertical  flame  spread  rates.  Since  the  ratio  of  length  to  solid  phase  time 

2/3 

(velocity)  should  increase  as  p  in  the  modeling  scheme,  flame  spread  rates 

2/3 

are  corrected  ftt  pressure  by  the  p  factor  shown  in  Figures  8-13.  As 
expected,  correlation  of  rates  of  change  of  flame  heights  are  far  less  satis¬ 
factory  than  the  original  flame  height  correlations.  Only  the  PMMA  corner 
configuration  still  yields  favorable  modeling  success  when  upward  spread  rates 
are  examined. 

Figures  8  and  9  show  clearly  that  modeling  of  upward  flame  spread  rate  on  a 
PMMA  wall  improves  as  air  pressure  increases  from  11.2  to  31.6  atmosphere. 

This  phenomenon,  which  has  been  discussed  previously,  is  due  to  the  competing 
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FIGURE  7  CORRELATION  OF  FLAME  HEIGHT  ON  PINE -WOOD  WALL  CORNERS 
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FIGURE  8  CORRELATION  OF  VERTICAL  FLAME  SPREAD  RATE  UP  PMMA  WALLS 
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effects  of  low  solid  surface  radiative  heat  loss  and  flame  radiative  satura¬ 
tion.  The  former  effect  is  more  important  at  11.2  atmoshperes,  leading  to 
spread  rates  higher  than  expected  from  the  modeling  scheme  while  flame  satura¬ 
tion  leads  to  reduced  spread  rates,  closer  to  the  expected  values  at  31.6 
atmospheres . 

Upward  flame  spread  rates  on  the  PMMA-inert  laminate  wall  are  correlated  in 
Figure  10.  Here,  the  approach  to  thermally  thick  behavior  as  elevated  pres¬ 
sure  increases  can  be  seen  more  clearly  than  was  the  case  for  the  flame  height 
data  in  Figure  4. 

For  the  wall-corner  configurations,  modeling  of  upward  flame  spread  rate  is 
only  accurate  for  the  PMMA  fuel,  as  is  evident  in  Figures  11-13.  Spread  rates 
for  particle-board  are  less  than  expected  at  11.2  and  20.5  atmospheres  but 
modeling  for  pine-wood  is  reasonably  good  at  these  pressures  (see  Figures  12 
and  13).  At  31.6  atmospheres,  scaled  spread  rates  for  both  particle-board  and 
pine-wood  are  much  less  than  corresponding  one-atmosphere  values. 

2.3.3  FLAME  WIDTH.  Data  on  the  maximum  width  of  the  flame  zone  (or  pyrolysis 
zone  for  the  case  of  the  PMMA  wall)  as  a  function  of  time  are  correlated  in 
Figures  14-19.  For  the  PMMA  wall  (see  Figure  14),  the  data  correlation  is 
quite  good  until  the  flame  extends  above  the  top  edge  of  the  model  and  full- 
scale  fuel  walls  at  t  *  600  and  800  seconds,  respectively.  This  degree  of 
modeling  success  is  somewhat  better  than  previous  pressure  modeling 
results^ perhaps  because  of  the  lack  of  side-walls  in  the  present  fuel 
configuration.  There  still  seems  to  be  a  strong  tendency  for  the  scaled, 
lateral  fire  growth  on  the  model  walls  to  be  more  rapid  than  that  at  one- 
atmosphere. 

This  same  tendency  is  evident  for  the  PMMA-PMMA  laminate,  data  for  which  are 
correlated  in  Figure  15.  The  lateral  spread  rate  for  both  the  model  and  full- 
scale  laminates  are  seen  to  be  somewhat  greater  than  that  for  a  uniform  PMMA 
wall.  With  the  PMMA-inert  laminate,  the  flame  widths  shown  in  Figure  16 
clearly  approach  those  characteristics  of  a  uniform  PMMA  wall  as  ambient 
pressure  increases  and  the  face  layer  becomes  thermally  thick.  As  a  result, 
scaled  lateral  spread  rates  for  the  model  laminates  are  less  than  the  values 
measured  for  the  full-scale  PMMA-inert  laminate. 
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VERT.  FLAME  SPREAD  RATE,  Vfy/p)  3  Cmm/e] 


FIGURE  13  CORRELATION  OF  VERTICAL  FLAME  SPREAD  RATE  UP  PINE-HOOD 
HALL  CORNERS 
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Figures  17-19  contain  flame  width  data  for  the  corner  fire  as  seen  by  an 
observer  on  the  bisector  of  the  90°  corner  angle.  Flame  width,  W,  is  there¬ 
fore  the  value  projected  onto  a  plane  normal  to  this  angle  bisector,  or  /  2 
times  the  horizontal  distance  on  the  fuel  surface  between  the  corner  apex  and 
each  side  of  a  symmetric  flame  front. 

Modeling  of  lateral  flame  spread  with  the  corner  configurations  is  seen  in 
Figures  17-19  to  be  reasonably  good.  In  fact,  the  correlation  of  data  is 
excellent  for  PMMA  fuel  (see  Figure  17).  Lateral  flame  spread  for  particle¬ 
board  is  modeled  well,  with  the  exception  of  the  31.6  atmosphere  test  The 
smaller  lateral  flame  widths  in  this  case  (see  Figure  18),  while  consistent 
with  the  mass  loss  data  shown  in  Figure  24,  may  simply  represent  data 
scatter.  There  appears  to  be  a  similar  type  of  behavior  for  the  pine-wood 
data  shown  in  Figure  19. 

2.3.4  MASS  LOSS  RATE.  Transient  measurements  of  fuel  mass  loss  during  upward 
and  lateral  fire  spread  are  fit  with  a  cubic  polynomial  least-6quares  regress¬ 
ion.  Differentiation  of  this  polynomial  fit  yields  the  fuel  mass  loss  rate, 
m  .which  equals  the  rates  of  soot  production  plus  fuel  gasification.  Since 
mass  flux  (kg/m  .s)  at  homologous  locations  should  increase  as  p  and  burn¬ 
ing  areas  decrease  as  p ~1*^  in  the  pressure  modeling  scheme,  total  mass  loss 

-2/3 

rate  should  decrease  as  p  .  Mass  loss  rates  are  therefore  corrected  for 

2/3 

pressure  by  use  of  a  p  factor  in  Figures  20-25  in  order  to  correlate  the 
data. 

Results  for  the  PMMA  wall  fire  are  shown  in  Figure  20.  Measurements  with  the 
smallest  model  fuel  at  31.6  atm  pressure  are  not  available  because  too  high  a 
load  system  sensitivity  is  required.  Up  to  the  time  of  extinguishment  of  the 
full-scale  wall,  modeling  of  mass  loss  rate  is  reasonable.  Much  of  the 
discrepancy  in  the  correlation  is  probably  due  to  the  lack  of  modeling  ofthe 
lateral  spread  process  after  a  test  time  of  800-900  seconds  (see  Figure  14). 

In  Figure  21 ,  the  mass  loss  rate  of  the  full-scale  PMMA-PMMA  laminate  is  seen 
to  be  somewhat  higher  than  that  of  a  uniform  PMilA  wall  of  the  same  size, 
probably  because  of  a  slightly  greater  lateral  flame  spread  (see  Figure  15)  in 
the  former  case.  The  greater  scaled  flame  heights  and  widths  on  the  model 
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FIGURE  19  CORRELATION  OF  MAXIMUM  PROJECTED  FLAME  WIDTH  ON 
PINE -WOOD  CORNERS 
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FIGURE  20  CORRELATION  OF  MASS  LOSS  RATE  FOR  PMMA  WALLS 
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laminate  fuels  (see  Figures  3  and  IS)  lead  to  the  higher  scaled  mass  loss 
rates  at  elevated  pressure  for  the  laminate.  Although  the  lack  of  data  corre¬ 
lation  at  later  test  times  is  due  to  he  fact  that  oversized  model  fuels  are 
used,  the  early  divergence  of  the  model  dat.  flame  dimensions  (and  hence 
mass  loss  rate)  is  difficult  to  explain. 

Equally  difficult  to  explain  is  the  high  scaled  mass  loss  rate  associated  with 
the  model  PMMA-Inert  laminate  at  20. 5  atm,  as  shown  in  Figure  22.  The  behav¬ 
ior  of  the  model  data  at  11.2  atm  and  at  31.6  atm  is  quite  reasonable  since 
the  model  laminates  are  approaching  the  thermally  thick  condition  of  the 
uniform  PMMA  wall. 

The  mass  loss  rate  correlations  for  both  the  PMMA  and  particle-board  wall- 
corners,  shown  in  Figures  23  and  24,  reflect  the  flame  spread  modeling  success 
for  these  configurations.  However,  such  is  not  the  case  for  the  pine-wood 
fuel,  as  seen  in  Figure  25.  It  is  apparent  that  the  low  mass  loss  rates  at 
one-atmosphere  must  be  due  to  flame  extinction  just  after  fire  spread.  Since 
char  formation  at  elevated  pressures  does  not  lead  to  a  similar  extinction 
phenomenon,  much  higher,  scaled  burning  rates  are  attained. 

Ill 

ANALYSIS 

3.1  LAMINATED  FUELS 

3.1.1  BACKGROUND .  Rigorous  application  of  the  pressure  modeling  technique 
would  require  that  the  dimensions  of  all  layers  of  a  laminated  fuel  should  be 
decreased  as  the  -2/3  power  of  absolute  air  pressure.  This  is  usually  not 
practical.  However,  many  real  materials  (perhaps  most)  are  composed  of  a  thin 
face  (or  wearing)  layer  laminated  to  much  thicker,  "backing"  layers.  During 
fire  spread  on  such  materials,  gasification  may  be  confined  to  the  face  layer 
alone  or  regression  through  the  face  layer  to  a  backing  layer  may  occur. 
Involvement  of  more  than  one  backing  layer  in  gasification  is  unlikely  during 
upward  flame  spread  on  vertical  walls  2-3  m  high  at  one  atmosphere.  In  fact, 
the  entire  thermal  wave  could  be  confined  to  the  thin  face  layer  during  upward 
spread  at  one  atmosphere  if  an  exposure  fire  imposed  a  sufficiently  large 

external  flux  on  the  wall  material. 


MASS  LOSS  R 


FIGURE  22  CORRELATION  OF  MASS  LOSS  RATE  FOR  PhMA- INERT  LAMINATE  WALLS 
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FIGURE  24  CORRELATION  OF  MASS  LOSS  RATE  FOR  PARTICLE -BOARD  WALL 
CORNERS 
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FIGURE  25  CORRELATION  OF  MASS  LOSS  RATE  FOR  PINE-WOOD  WALL  CORNERS 
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3.1.2  BEHAVIOR  OF  THERMAL  WAVE  DURING  FIRE  SPREAD.  Assume  fuel  gasification 
is  confined  to  the  "face"  layer  of  a  laminated  material.  By  definition  the 
face  layer  would  be  thermally  thick  if  the  entire  thermal  wave  were  also 
confined  to  this  layer  during  upward  spread.  Pressure  modeling  of  such  a 
"thermally  thick"  layer  would  require  no  modification  of  layer  thickness. 

A  much  more  likely  scenario  is  that  the  thermal  wave  extends  into  the  backing 

layer  during  one  atmosphere,  upward  spread.  Since  the  thermal  wave  thickness^7^ 

is  of  the  order,  a/V,  where  a  is  the  thermal  diffusivity  (X/pc)  and  V  is  the 

fuel  regression  rate,  the  increase  in  fuel  regression  rate  as  the  2/3  power  of 

absolute  air  pressure  should  lead  to  a  decrease  in  thermal  wave  thickness  as 
-2/3 

p  .  For  sufficiently  high  air  pressure,  p,  the  entire  thermal  wave  could 
be  confined  to  the  face  layer  of  the  model  fuel  if  the  face  layer  thickness 
were  not  reduced  from  the  full-scale  value.  A  face  layer  of  PMMA  about  2.7  mm 
thick,  for  example,  would  probably  contain  the  entire  thermal  wave  at  30 
atmospheres  during  upward  spread  on  a  wall  0.25  m  high  (modeling  a  2.5  m  i igh 
full-scale  wall). 

The  implications  of  this  decrease  in  thermal  wave  coverage  can  be  determined 
from  a  simplified  one-dimensional  analysis  of  flame  spread  rates.  Assume  that 
the  thermal  wave  extends  over  a  face  layer  (subscript  "1")  of  thickness,  d, 
and  penetrates  to  a  distance,  6  ,  from  the  face  layer  into  the  backing  layer 
(subscript  "2").  Furthermore,  assume  that  flame  spread  occurs  when  a  net 
flame  heat  flux,  q"  ,  raises  the  surface  temperature  of  the  face  layer  from 
ambient,  ,  up  to  the  pyrolysis  temperature,  T  during  a  time  interval. 

At  .  If  the  temperature  of  the  heated  backing,  is  T  at  the  face  layer 
boundary  and  a  distance,  6  ,  from  the  face  layer  when  pyrolysis  first 
occurs,  then  from  energy  conservation: 

* 

T  +  T  * 

q"  At  -  PjCjd  (-*-2 - Too>  +  P2C26  (T  "  T»)/2  (1) 

Because  of  thermal  conduction  across  the  face  layer: 

q"  *  (T  -  T* )/d 

Across  the  heated  portion  of  the  backing,  thermal  conduction  gives: 
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(2) 


(3) 


q"  -  a2  (T  -  TJ/6 


From  Eqs.  (2)  and  (3): 


T  -T  =  T  -T  -  q”  d/A, 

«  p  ®  1 

A„ 

«  “  ~  IT  "  T»  ~  d/V 

... 


(4) 

(5) 


Substitution  of  Eqs.  (4)  and  (5)  into  Eq.  (1)  yields  an  expression  for  the 
time  interval.  At  ,  required  to  heat  the  fuel  surface  to  the  pyrolysis  temper¬ 
ature  : 


p.c.d  (T  -  T  -  q"  d/2A  )  p,c0A,  . 

- E - = - L_  +  _2_2J_  [T  _T  _  ^  d/x  ,2 

....  2  p  “>  1 


2q" 


(6) 


The  flame  spread  race  is  inversely  proportional  to  this  heating  time 
inverval.  For  upward  fire  spread,  it  is  known^®^  that  Vp  is  also  directly 
proportional  to  the  current  pyrolysis  height,  xp.  Thus, 

V  “x  /At  (7) 

P  P 

When  the  fuel  consists  of  a  uniform,  thermally  thick  material  (no  face  layer), 
then  d  *  0  in  Eq.  (6).  The  spread  velocity,  from  Eqs.  (6)  and  (7),  becomes: 

V  «  2  x  q"2/pcX  (T  -  T  )2  (8) 

p  p  P 

When  the  fuel  consists  of  a  thermally  thin  face  layer  alone,  with  no  backing 
material,  then  Tp  *  T*  in  Equation  4  and  p^c^A^  “  ®  *n  ^9*  The  spread 

velocity,  from  Eqs.  (6)  and  (7),  is  then  given  by: 

V  <*  x  q"/pcd  (T  -  T  )  (9) 

P  P  P  ® 

With  a  laminated  fuel  consisting  of  both  a  face  layer  and  backing,  flame 
spread  velocity  is  derived  from  the  complete  expression  in  Eq.  (6)  as 
follows : 
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PlCld  ^Tp  ~  T®  ~  ^  d/ 2 A  ^  >  p2c2X2 


-  +  T^T  (TP '  T«  - «"  d/V  J 

^9  nr.) 


The  preceding  relations  can  be  used  to  show  qualitatively  the  relative  magni¬ 
tudes  of  flame  spread  rates  for  various  types  of  materials.  For  Instance,  the 
PMMA  fuel  used  for  the  laminated  wall  fires  in  the  present  study  has  a  pyroly¬ 
sis  (or  vaporization)  temperature  at  one-atmosphere  close  to  636  K  (see 
Ref.  19)  and  other  thermal  properties  given  in  Table  1.  With  an  assumed  net 
heat  flux  to  the  PMMA  of  104  W/m  (1  W/cm  the  ratio  of  upward  spread  rate 
to  the  pyrolysis  height,  Vp/x^,  is  3.6  times  greater  for  an  isolated  1  mm 
thick  layer  (Eq.  (9))  than  for  a  thermally  thick  slab  (Eq.  (8)).  Use  of 
Eq.  (10)  shows  that  the  same  1  mm  PMMA  layer  backed  by  the  inert  ceramic 
material  should  support  an  upward  spread  rate  nearly  as  large  as  the  rate  with 
the  unbacked  layer  alone,  or  a  ratio  3.55  times  greater  than  for  the 

thermally  thick  slab. 

An  examination  of  Figures  4  and  10  allows  the  measured,  full-scale  upward 
spread  rates  on  thermally  thick  PMMA  and  the  PMMA-Inert  laminate  to  be  compar¬ 
ed  for  the  same  values  of  flame  height.  For  instance,  at  a  flame  height,  x^ , 
of  0.6  m  (t  *  700g  for  the  uniform  slab,  t  =  350g  for  the  laminate),  the 
upward  spread  rate  on  the  laminate  is  about  2.2  times  greater  than  that  on  the 
uniform  wall.  Assuming  that  the  0.6  m  flame  height  depends  only  on  pyrolysis 
height  for  the  PMMA  fuel  walls,  the  ratio,  Vp/Xp,  for  the  laminate  should  also 
be  about  2.2  times  larger  than  that  for  the  uniform  wall.  Use  of  the  simpli¬ 
fied  analysis  with  the  actual  PMMA  face  layer  thickness,  d,  of  3  mn  on  the 
inert  backing  yields  a  ratio,  Vp/Xp,  which  is  about  1.5  times  higher  than  the 
value  corresponding  to  thermally  thick  PMMA.  The  analysis  thus  seems  to  give 
a  good  qualitative  description  of  the  dependence  of  spread  rates  on  fuel 
thermal  characteristics. 


Solid  phase  thermal  conduction  and  upward  spread  rate  for  laminated  materials 
at  elevated  ambient  pressures  can  also  be  studied  with  the  preceding  rela¬ 
tions.  From  Eq.  (5),  it  is  seen  that  the  depth,  6  ,  of  thermal  wave  penetra¬ 
tion  into  the  backing  material  decreases  sharply  as  the  net  heat  flux  to  the 
fuel  increases  with  pressure  (roughly  as  p  '  ).  If  the  face  layer  thickness. 


d,  is  reduced  as  p-'^,  q"  d  will  remain  constant  and  then,  from  Eq.  (5),  6 
will  simply  decrease  as  the  required  p-^^. 

A  constant  face  layer  thickness,  however,  is  shown  by  Eq.  (5)  to  result  in  a 

6  of  zero  when  q"  »  X,  (T  -  T  )/d  .  For  a  3  m  thick  face  layer  of  PMMA, 

1  P  4  2  ? 

this  condition  is  satisfied  when  q"  -  2.35  x  10  W  a  (or  2.35  W/cnr),  at 

which  point  there  would  be  no  penetration  of  the  thermal  wave  into  the  backing 

material.  Such  an  increase  in  heat  flux  might  occur  at  a  pressure  of  only 

A  o 

3.6  atmospheres  if  the  net  flux  at  one-atmosphere  is  1(T  W/m  .  At  higher 

ambient  pressures,  the  constant  thickness  face  layer  then  becomes  thermally 

2 

chick  with  a  spread  rate  given  by  Eq.  (8)  and  V  /x  increasing  as  q"  ,  or 

For  a  face  layer  thickness  varying  as  it  can  be  shown  that  the 

spread  rate  relation  in  Eq.  (10)  also  yields  a  value  of  vp/xp  increasing  as 

and  (as  noted  previously)  always  substantially  greater  than  that  given  by 

Eq.  (8)  if  the  backing  is  the  inert  ceramic.  Thus,  lower  scaled  spread  rate 

parameters,  V  /x  ,  will  be  measured  for  the  constant  thickness  model  than  for 
?  P 

the  full-scale  laminate. 

3.1.3  PRESSURE  MODELING  LAMINATED  FUELS  A  simplified  heat  transfer  analysis 
shows  that  scaled  spread  rates  much  lower  than  full-scale  values  would  be 
predicted  by  pressure  modeling  if  1)  the  pel  of  the  constant  thickness  face 
layer  of  a  laminate  is  much  greater  than  that  of  the  backing  layer  and  2)  the 
face  layer  thickness  is  sufficiently  small  yet  still  contains  most  of  the 
thermal  wave  at  elevated  pressures.  Such  a  situation  could  probably  be 
expected  for  most  real  composite  materials  since  the  face  layer  is  usually 
quite  thin  and  dense  (for  wear  resistance)  while  the  backing  may  be  a  low 
density  foam  or  honeycomb. 

This  modeling  behavior  could  be  corrected  by  reducing  the  face  layer  thickness 
as  p-^^.  Such  reductions  may  not  always  require  fabrication  of  new  laminated 
materials  but  possibly  a  machining  operation  (grinding  or  sanding)  of  the 
surface  of  the  small-scale  fuel. 

Assume  that  gasification  extends  into  the  backing  layer  during  upward  fire 
spread  at  one  atmosphere.  In  this  case,  accurate  modeling  would  require  a 
reduction  (by  machining)  of  surface  layer  thickness  so  that  the  model  contains 
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the  correct  amount  of  surface  layer  fuel.  Accurate  machining  would  be  diffi¬ 
cult,  however.  Alternatively,  two  cases  could  be  examined:  model  spread 
rates  both  with  and  without  the  entire,  unmodified  face  layer  can  be  obtained 
and  the  more  hazardous  result  used  to  characterize  the  material. 

3.2  CHARRING  FUELS 

The  heat  and  mass  transfer  processes  within  a  charring  fuel  can  be  modeled 
numerically  to  determine  how  pyrolysis  reactions  will  influence  pressure 

modeling  success.  Such  calculations  have  been  performed  using  the  SPYVAP 
(9) 

computer  code  for  transient,  one-dimensional  thermal  conduction  and  pyroly¬ 
sis  with  one-step  Arrhenius  kinetics.  This  numerical  procedure  is  documented 
in  detail  in  Appendix  A,  which  is  taken  from  Reference  9.  Parameters  for  the 
numerical  procedure  are  listed  in  Table  III  and  computed  results  are  given  in 
Figures  26-33.  An  external  radiant  flux  is  assumed  to  have  a  constant  value 
beginning  at  the  start  (t  **  0)  of  the  transient  pyrolysis  process  to  simulate 
the  presence  of  a  flame.  Convective  heat  flux  also  is  allowed  due  to  an 
assumed,  constant  flame  temperature  of  1350  K.  Unless  noted  otherwise,  the 
optical  depth  (kLm),  material  thickness  and  heat  transfer  coefficients  are  all 
scaled  according  to  the  pressure  modeling  scheme.  Values  of  optical  depth  at 
one-atmosphere  (khm)0  are  selected  to  yield  a  flame  radiant  (exposure)  flux  of 
either  20  or  40  kW/m^  with  a  1350  K  flame  temperature,  as  noted  on  each  of 
Figures  26-',3.  At  elevated  pressures,  the  radiant  exposure  flux,  q^  ,  is 
given  by  the  following  expression  (see  Reference  1): 

q"  -  a  (1350  K) ^  [1  -  exp  (-  (^ — )2^  (kL  )  )] 
r  p  m  o 

ro 

where  a  is  the  Stefan-Boltzmann  constant. 

As  shown  in  Table  III,  the  char-forming  wood  fuels  are  assumed  to  have  speci¬ 
fic  heats  and  thermal  conductivities  which  are  functions  of  both  the  local, 
solid-phase  temperature,  T [ K ] ,  and  the  relative  amounts  of  char  and  virgin 
material  at  each  instant.  These  relations  and  the  remaining  kinetics  parame¬ 
ters  were  obtained  from  studies  by  Kung^®’^\  from  previous  SPYVAP  calcula- 

(12) 

tions  performed  by  Tamanini  and  from  extrapolations  of  measurements 
reported  in  References  13-18. 
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TABLE  III 


3.2.1  FUEL  MASS  FLUX  In  Figures  26-29,  computed  values  of  fuel  mass  flux  at 
the  "front"  face  of  the  fuel  (see  Table  III)  are  correlated  for  one-atmosphere 
and  two  elevated  pressure  conditions.  Mass  flux  is  corrected  for  pressure  as 
shown  in  these  figures  because  m"  should  increase  as  p  '  J  in  the  modeling 

scheme.  As  noted  before,  fuel  thickness  is  generally  reduced  for  the  calcula- 

-2/3 

tions  as  p  from  the  full-scale  values  given  in  Table  III.  Although  such  a 
reduction  in  thickness  is  not  made  for  the  char-forming  fuels  in  the  actual 
elevated  pressure  experiments,  the  numerical  solution  technique  requires  the 
reduction  in  order  to  have  an  adequate  number  of  grid  points  within  the 
thermal  wave. 

Mass  flux  from  a  pinewood  fuel,  simulated  with  the  parameters  given  in 
Table  III,  is  shown  in  Figures  26-28  for  an  exposure  radiant  flux  at  one- 
atmosphere  of  20  or  40  kW/m  .  Initially,  there  is  a  rapid  rise  in  fuel  mass 
flux  to  a  peak  value  due  to  the  radiant  exposure.  This  fuel  “pulse"  is  then 
followed  by  a  decay  period  because  of  the  increasing  thickness  of  the  char, 
which  insulates  the  virgin  fuel. 

It  appears  that  the  higher,  most  widely  accepted  value^*^  of  Arrhenius 
activation  energy  for  fuel  pyrolysis  yields  the  best  pressure  modeling  of  the 
pinewood  fuel  mass  flux.  Modeling  errors  also  are  reduced  for  the  higher 
exposure  flux.  The  actual  radiant  flux  to  the  charring  fuels  in  the  corner 

O 

configuration  is  likely  to  be  close  to  this  40  kW/m*  value  due  to  radiant 
emission  from  the  adjacent,  high  temperature  fuel  surfaces  (see  Figures  30- 
33).  In  any  case,  it  is  probable  that  adequate  pressure  modeling  of  fuel  mass 
flux  would  lead  to  corresponding  success  in  modeling  the  actual  flame  heat 
transfer,  and  hence  the  fire  spread  process. 

In  Figure  29,  the  calculated  mass  flux  corresponding  to  a  higher  density, 
char-forming  fuel,  such  as  particle-board,  is  correlated.  The  predicted 
success  of  pressure  modeling  in  this  case  is  not  as  good  as  that  for  the  lower 
density  pinewood.  This  prediction  seems  to  be  sustained  by  the  correlation  of 
flame  height  measurements  in  Figures  6  and  7. 
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FIGURE  27  CALCULATION  OF  ONE -DIMENSIONAL  PINE-WOOD  PYROLYSIS:  SCALED 
VAPOR  MASS  FLUX  DUE  TO  20  kW/tn2  SCALED  HEAT  FLUX  WITH  LOW 
ACTIVATION  ENERGY  IN  PYROLYSIS  KINETICS 
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FIGURE  28  CALCULATION  OF  ONE -DIMENSIONAL  PINE-WOOD  PYROLYSIS:  SCALED 

VAPOR  MASS  FLUX  DUE  TO  40  kW/m2  SCALED  HEAT  FLUX  WITH  STANDARD 
THERMAL  AND  KINETICS  PROPERTIES 
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As  described  in  Section  2.2,  the  cellulosic  fuels  at  one-atmosphere  undergo  a 
flame  extinction  process  after  fire  spread  proceeds  some  distance  from  the 
corner  apex.  Empirical  studies  with  a  vareity  of  materials  have  shown  that 
such  flame  extinction  occurs  when  fuel  mass  flux  drops  below  the  2  to  4  g/m  .s 
range.  The  calculated  mass  flux  at  one-atmosphere  is  seen  in  Figures  26  and 
29  to  fall  below  the  3  g/m  .s  level  some  380  and  560  seconds  after  initial 
radiant  exposure  of  pinewood  and  particle-board,  respectively.  Coincidental¬ 
ly,  the  measured  rates  of  increase  in  flame  height  on  pinewood  and  particle¬ 
board  corners,  shown  in  Figures  13  and  12  respectively,  drop  to  zero  at  about 
these  same  times.  Lateral  flame  spread  on  the  full-scale  fuels  stops  as  a 
direct  result  of  the  extinction  phenomenon  at  respective  times  after  ignition 
of  360  and  690  seconds  for  pinewood  and  particle-board,  as  shown  in  Figures  19 
and  18. 


3.2.2  FUEL  SURFACE  TEMPERATURE  The  calculated  temperature  of  the  "front" 
face  of  the  fuel  during  radiant  exposure  and  fuel  pyrolysis  is  shown  in 
Figures  30-33.  In  all  cases,  the  fuel  surface  temperature  reaches  a  steady 
value  within  about  100  seconds  after  the  start  of  radiant  exposure.  Predicted 
values  of  char  surface  temperature  at  one-atmosphere  are  about  800  K  (900  K 
with  a  A0  kW/m2  flux)  for  the  cellulosic  fuels. 

At  elevated  pressures,  the  computed  surface  temperatures  for  the  char-forming, 
cellulosic  fuels  are  quite  high,  up  to  1200  K  at  31.6  atm.  Surface  tempera¬ 
ture  itself  is  thus  not  being  pressure  modeled  (i.e.  preserved),  but  the 
resultant  enhancement  of  radiant  heat  loss  from  the  surface  does  allow  for 
pressure  modeling  of  the  net  heat  flux  to  the  fuel.  If  both  heat  loss  due  to 
surface  reradiation  and  heat  gain  due  to  flame  heat  flux  increase  roughly  as 
required  by  the  modeling  scheme,  then  the  net  radiant  heat  flux  and  hence  the 
fuel  mass  flux,  should  be  modeled  well  (assuming  convective  heat  transfer  is 
relatively  unimportant).  The  predicted  heat  loss  due  to  surface  emission  is 
in  fact  seen  to  increase  somewhat  less  than  p*'  (a  factor  of  3.6  rather  than 
5)  between  one  and  11.2  atm  for  pinewood  and  particle-board  fuels.  This  may 
be  the  reason  for  the  success  in  modeling  flame  growth  on  the  cellulosic  wall- 
corners. 
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FIGURE  30  CALCULATION  OF  ONE-DIMENSIONAL  PINE-WOOD  PYROLYSIS;  FUEL 
SURFACE  TEMPERATURE  DUE  TO  20  kW/m2  SCALED  HEAT  FLUX  WITH 
STANDARD  THERMAL  AND  KINETICS  PROPERTIES 
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CALCULATION  OF  ONE -DIMENSIONAL  PINE-WOOD  PYROLYSIS:  FUEL 
SURFACE  TEMPERATURE  DUE  TO  20  kW/ra2  SCALED  HEAT  FLUX  WITH 
LOW  ACTIVATION  ENERGY  IN  PYROLYSIS  KINETICS 
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FIGURE  32  CALCULATION  OF  ONE -DIMENSIONAL  PINE-WOOD  PYROLYSIS:  FUEL 
SURFACE  TEMPERATURE  DUE  TO  40  kW/m2  SCALED  HEAT  FLUX  WITH 
STANDARD  THERMAL  AND  KINETICS  PROPERTIES 
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FIGURE  33  CALCULATION  OF  ONE-DIMENSIONAL  PARTir-  p  » 
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3.2.3  CHAR  PRODUCTION  Another  important  output  of  the  numerical  calculations 
is  the  predicted  thickness  of  the  char  layer  which  develops  at  a  given  loca¬ 
tion  after  fuel  pyrolysis  is  complete.  When  pine-wood  fuel  is  simulated  (high 
activation  energy)  with  a  40  kW/m  exposure,  this  char  layer  grows  to  about 
43%  of  the  total  fuel  thickness  after  800  seconds  of  radiant  exposure  at  one- 
atmosphere.  Corresponding  char  fractions  of  scaled  fuel  thickness  at  11.2  and 
31.6  atmospheres  are  62%  and  71%  for  the  same  scaled  exposure  times.  It 
appears  from  these  data  and  from  results  with  a  20  kW/m  exposure  that  the 
scaled  char  thickness  increases  over  the  value  at  one-atmosphere  by  a  factor 
of  (p/pQ)®*^  when  che  high  activation  energy  kinetics  is  used.  While  any 
such  pressure  dependence  represents  a  departure  from  the  modeling  scheme,  the 
small  increase  in  relative  char  thickness  helps  to  reduce  fuel  mass  flux,  and 
thus  enhances  the  effect  of  the  calculated  increase  in  fuel  surface  tempera¬ 
ture  at  elevated  pressure.  Together,  the  insulating  property  of  the  excess 
char  formation  and  the  heat  loss  from  a  high  temperature  surface  act  to  keep 
fuel  mass  flux  increasing  at  close  to  the  required  2/3  power  of  pressure. 

It  should  be  noted  that  use  of  the  low  activation  energy  kinetics  leads  to  a 
calculated  char  fraction  (of  fuel  thickness)  which  is  nearly  independent  of 
pressure  (about  33%  at  20  kW/m^). 


IV 

SUMMARY  OF  RESULTS 

Measurements  of  time-resolved  flame  dimensions,  upward  spread  rates  and  fuel 
mass  loss  rates  for  full-scale  and  model  configurations  burning  at  ambient  and 
elevated  air  pressures,  respectively,  yield  the  following  results: 

1.  Pressure  modeling  is  sufficiently  accurate  for  the  prediction  of  fire 
growth  from  a  point  ignition  on  a  uniform  PMMA  wall  when  both  upward  and 
lateral  flame  spread  processes  occur. 

2.  The  behavior  of  the  flame  spread  process  at  elevated  air  pressures  for 
walls  composed  of  a  face  layer  of  PMMA  with  a  thick  nonburning  backing  layer 
is  not  completely  consistent  with  a  simplified  analysis  oi  thermal  conduction 
processes.  When  the  backing  layer  is  also  PMMA,  spread  rates  at  elevated 
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pressures  are  much  greater  than  expected,  but  spread  rates  are  more  consistent 
with  heat  transfer  theory  when  the  backing  is  a  low  density  ceramic. 

3.  Pressure  modeling  of  fire  growth  in  a  wall  corner  configuration  is  quite 
accurate  for  cellulosic,  char-forming  materials  and  for  PMMA.  The  cellulosics 
at  one  atmosphere  exhibit  a  flame  extinction  phenomenon  after  significant 
lateral  flame  spread  that  is  not  observed  at  elevated  pressures. 

4.  Thermally  thick,  rigid,  high  density  polyurethane  foam  in  a  corner  config¬ 
uration  will  not  support  significant  flame  spread  at  one-atmosphere  but  will 
at  elevated  pressures  with  a  properly  scaled,  small  PMMA  ignition  source. 

This  behavior  is  perhaps  due  to  excessive  radiant  heat  loss  from  the  char  and 
the  intumescent  character  of  the  char  at  one-atmosphere.  Gas  phase  chemical 
kinetics,  which  may  be  the  most  important  factor  in  the  initiation  of  flame 
spread  on  the  full-scale  foam,  is  clearly  not  modeled. 

5.  A  simplified  analysis  of  thermal  conduction  in  a  laminated  material  is 
used  to  show  how  flame  spread  rates  are  affected  by  the  thermal  properties  of 
the  face  layer  and  backing  at  both  one-atmosphere  and  at  elevated  pressures. 

6.  A  numerical  technique  is  used  to  predict  one-dimensional,  transient  fuel 
mass  flux,  fuel  surface  temperature  and  char  thickness  during  material  expo¬ 
sure  to  a  prescribed  radiant  (and  convective)  heat  flux.  Calculated  results 
show  that  reasonable  pressure  modeling  of  flame  spread  rates  should  be  expect¬ 
ed  for  the  cellulosic  fuels  due  to  increases  in  surface  temperature  and  char 
production  for  conditions  simulating  elevated  air  pressure.  Predicted  times 
for  extinction  of  these  fuels  at  one-atmosphere,  due  to  char  buildup,  are  in 
agreement  with  the  observed  times. 


V 

CONCLUSIONS 

1.  Pressure  modeling  of  three-dimensional  fire  spread  on  uniform  walls  and 
wall  corners  has  now  been  validated  for  PMMA  and  for  wood  fuels.  It  has  not 
yet  been  established  that  the  modeling  technique  is  also  valid  for  predicting 
fire  growth  on  other  charring  fuels  in  corner  configurations.  Results  from 
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this  and  previous  studies  have  shown  that  in  general,  the  complex  process  by 
which  self-sustained  fire  spread  is  initiated  is  not  pressure  modeled.  Such 
fire  spread  initiation  occurs  readily  at  elevated  pressures  because  surface 
radiant  heat  loss  and  the  action  of  gas  phase  chemical  retardents  cannot  be 
modeled.  With  wood  and  PMMA  fuels,  fire  spread  rates  on  wall-ceiling  corners 
should  also  be  predictable  by  pressure  modeling,  based  on  previous  work^®^ 
with  ceiling  channel  configurations. 

2.  It  appears  that  much  more  experimental  information  is  needed  before 
pressure  modeling  can  be  used  to  predict  fire  growth  on  practical  composite 
materials.  At  present  the  thickness  of  all  components  (including  adhesives) 
within  the  thermal  wave  developed  during  fire  spread  must  be  modified  accord¬ 
ing  to  the  pressure  modeling  scheme.  However,  radiant  exposure  in  real 
enclosure  fires  may  well  be  sufficiently  high  (2-4  W/cm  )  to  confine  thermal 
wave  penetration  to  a  surface  layer  of  fuel  during  fire  spread.  Pressure 
modeling  of  such  a  fire  spread  process  would  then  be  accurate  without  any 
modification  of  the  laminated  material. 

3.  Pressure  modeling  can  serve  as  a  scientific  tool  for  studying  fire  growth 
on  idealized  charring  and  laminate  fuels  in  a  variety  of  configurations.  In 
this  way,  the  potential  for  fire  growth  as  a  function  of  fuel  geometry  can  be 
determined.  However,  further  study  is  needed  to  see  if  pressure  modeling 
correctly  predicts  the  relative  hazard  of  different  fuel  compositions. 
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APPENDIX  A 

A  NUMERICAL  MODEL  FOR  ONE-DIMENSIONAL  HEAT 
CONDUCTION  WITH  PYROLYSIS  IN  A  SLAB  OF  FINITE  THICKNESS 


F.  Taman ini 

Factory  Mutual  Research  Corporation 

A.  1  INTRODUCTION 

The  purpose  of  this  appendix  is  to  document  a  procedure  for  computing  the  pro¬ 
files  of  temperature,  density  and  mass  flux,  as  well  as  the  surface  energy  fluxes 
associated  with  them,  in  a  slab  of  finite  thickness  undergoing  pyrolysis.  No  new 
physics  is  introduced  to  make  the  model  more  realistic  than  similar  models  de¬ 
veloped  by  other  researchers.  In  particular,  the  fundamental  equations  are 

1  2 

those  proposed  by  Kung  and  later  used  by  this  writer  ;  with  respect  to  those 
versions  of  the  model,  however,  the  current  procedure  offers  greater  flexibility 
and  a  few  additional  options. 

The  features  of  the  model  are: 

1)  Heat  conduction  Is  calculated  by  allowing  for  variable  thermal  proper¬ 
ties.  The  thermal  conductivity  (k)  and  the  specific  heat  (c^)  are  assumed  to  be 
linear  functions  of  the  local  temperature. 

2)  Pyrolysis  follows  a  first  order  Arrhenius  reaction:  the  thermal  de¬ 
composition  transforms  active  material  into  constant,  pre-assigned  fractions  of 
char  and  volatiles.  Before  pyrolysis  is  completed  the  solid  matrix  consists  of 
a  mixture  of  char  and  unpyrolyzed  active  material,  whose  thermal  properties  are 
obtained  by  linear  interpolation  of  the  property  values  pertaining  to  the  two 
components. 


"Stung,  H.C.,  "A  Mathematical  Model  of  Wood  Pyrolysis,"  Combustion  and  Flame, 
18,  185-195  (1972) 

2 

Tamanini,  F.,  "A  Study  of  the  Extinguishment  of  Vertical  Wood  Slabs  in 
Self-sustained  Burning  by  Water  Spray  Application,"  Combustion  Science  and 
Technology,  14,  1,2,3,  p.  1  (1976)  and  "Everything  You  Always  Wanted  To  Know 
From  a  Thermocouple  (in  a  fire  test).  But  Were  Afraid  To  Ask,"  Society  of 
Fire  Protection  Engineers,  Technology  Report  75-2,  (1975) 
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3)  Thermal  decomposition  contributes  to  the  local  energy  balance  through 
a  volume  generation  of  heat.  The  heat  of  pyrolysis  associated  with  that  process 
is  assumed  to  be  constant  at  a  reference  temperature. 

4)  Accumulation  of  volatiles  within  the  solid  matrix  is  neglected.  All 
the  gaseous  products  of  Che  decomposition  process  are  assumed  to  escape  toward 
either  or  both  surfaces  as  they  are  generated. 

5)  Convective  heat  transfer  between  the  volatiles  and  the  solid  matrix  is 
taken  into  account  by  postulating  perfect  thermal  contact. 

6)  Boundary  conditions  at  the  two  bounding  surfaces  can  be  specified  in 
terms  of  temperature  or  heat  flux.  If  the  temperature  is  prescribed,  the  model 
yields  the  heat  flux  and  vice  versa.  To  allow  for  a  situation  often  encountered 
in  practice,  the  model  can  also  use  as  a  boundary  condition  a  temperature -time 
history  at  a  location  inside  the  slab. 

The  computer  program  illustrated  here  consists  of  a  MAIN  program  and  two 
subroutines:  SPYVAP  (Slab  Pyrolysis  with  Variable  Properties)  and  OUTPUT.  The 
following  sections  concentrate  on  the  description  of  subroutine  SPYVAP,  which 
contains  the  main  machinery  of  the  model.  Subroutine  OUTPUT  calculates  secondary 
quantities  of  interest,  such  as  heat  flux  components,  mean  slab  density,  etc., 
but  its  main  function  is  to  do  just  what  its  name  implies.  The  version  of  MAIN 
reported  here  is  the  one  used  to  perform  the  calculations  discussed  in  Section  V 
of  this  report.  The  function  of  MAIN  is  to  supervise  the  calls  of  the  two  sub¬ 
routines  as  well  as  to  initialize  the  array  containing  the  temperature/heat  flux 
profiles  used  as  boundary  conditions  and  effect  input/output  of  the  initial  con¬ 
ditions.  Users  of  the  procedure  should  not  need  to  modify  SPYVAP  but  will  have 
to  rewrite  MAIN  to  adapt  it  to  their  particular  application.  Changes  in  OUTPUT 
may  also  be  necessary  to  satisfy  personal  aesthetic  requirements  or  special 
output  needs. 
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A. 2  MODEL  EQUATIONS  AND  BOUNDARY  CONDITIONS 

The  model  finds  solutions  to  the  problem  of  unsteady  heat  conduction  in  one 
dimension.  The  governing  equation,  with  the  appropriate  terms  to  account  for 
convective  heat  transfer  between  the  solid  and  the  volatiles  and  for  energy  re¬ 
lease  in  the  pyrolysis  process,  reads: 


3(p  h  ) 
s  s 

3t 


+  —  (M  h  ) 
8  g 


(1) 


where :  t  time , 

x  distance  from  the  front  surface, 
p  density, 

k  thermal  conductivity, 
h  enthalpy , 

M  mass  flux  of  volatiles  (positive  in  the  negative-x  direction) , 

8 

Q  heat  of  pyrolysis  (positive  when  reaction  is  exothermic) . 

The  subscripts  s  and  g  refer  to  the  solid  matrix  and  the  volatiles  (pyrolysis 
gases)  respectively.  The.  enthalpy,  h,  is  defined  as: 


where : 


T 

h  -  /  c  (T)  dT 

T  P 
o 


VT> 


O  * 

c +  c  (T  -  T  ) 
p  p  o 


(2) 

(2') 


The  mass  flux  of  pyrolysis  gases,  M 


g’ 


is  calculated  from: 


3M  dp 

— g _ §L 

3x  3t 


(3) 


As  for  determining  the  direction  of  the  migration  of  the  pyrolysis  gases,  the 
model  offers  two  options? 
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1)  The  volatiles  flow  in  the  direction  of  decreasing  densities  of  the  solid 
matrix. 

2)  The  outflow  of  volatiles  is  laminar  and  there  is  no  net  pressure  dif¬ 
ference  between  the  two  faces  of  the  slab. 

Since  in  laminar  flow  the  pressure  drop  is  proportional  to  the  flow  mean  velocity, 
the  following  condition  is  used  to  implement  the  latter  of  the  two  options: 

£ 

/  M  (x)  dx  ■  0  (A) 

0  8 

where  £  is  the  thickness  of  the  slab.  To  arrive  at  eq  (A)  changes  in  the  density 
of  the  volatiles,  as  well  as  changes  in  the  porosity  of  the  solid  matrix,  have  been 
neglected. 

After  a  certain  amount  of  pyrolysis  has  occurred,  part  of  the  solid  matrix  is 

char,  the  rest  is  as  yet  unpyrolyzed  active  material.  It  is  assumed  that  the 

density  of  the  active  material  (p  )  and  that  of  the  solid  matrix  (p  )  are  linearly 

fl  s 

related: 


Pg(t,x) 


fi 


(t,x)  +  pf 


(5) 


where  p^  and  pf  are  the  initial  and  final  densities  respectively. 

The  density  of  the  char  can  be  calculated  from  values  for  p  and  p  as  : 

as 


p  (t,x)  -  p  (t ,x)  -  P  (t ,x)  (6) 

c  s  a 

At  the  beginning  p  •  p  ■  p.  and  p  “0;  after  complete  pyrolysis  p  *  p  *  Pf 
fl  6  X  C  C  6  I 

and  p  -  0. 

a 

The  rate  of  pyrolysis  is  determined  by  using  a  first-order  Arrhenius 


reaction: 
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where  a  is  the  pre-exponential  factor,  E  is  the  activation  energy  and  R  the 
P 

gas  constant. 

Possible  recondensation  within  the  solid  is  not  taken  into  account.  It  should 

be  noted  that  eq  (7)  is  written  by  some  authors  with  p  -  p  instead  of  p  in  the 

right  hand  side.  Since  the  two  quantities  are  proportional  to  each  other 

(cf  eq  (5)),  the  difference  is  conceptually  not  too  important.  However,  the 

adoption  of  p  instead  of  p  -  p.  amounts  to  introducing  the  factor  (1  -  pf/p . ) 

which  must  be  taken  into  account  before  some  of  the  values  for  a  available  in 

P 

the  literature  can  be  used  in  the  model. 

The  contributions  from  the  char  and  the  active  material  to  the  energy 
content  of  the  solid  is  expressed  as: 


P  h 
s  s 


p  h  + 
a  a 


P  h 
c  c 


(8) 


Equation  (8)  can  be  used  to  obtain  an  expression  for  the  mean  specific  heat  of 

the  solid  (c  )  as  a  function  of  those  of  the  active  material  (c  )  and  the  char 
ps  pa 

(c  )  and  the  local  density  (p  )  : 

pc  s 


c 

ps 


(8’) 


where  p&  and  p ^  can  be  written  as  functions  of  pg  using  eq  (5). 

For  the  thermal  conductivity  kg,  a  linear  variation  with  density  is  assumed 

between  the  value  of  the  active  material  (k  )  and  that  of  the  char  (k  )  : 

a  c 


k 

s 


(9) 


Again,  with  the  aid  of  eq  (5)  pg  can  be  substituted  in  the  above  relationship  to 
and  p^.  Furthermore,  k&  and  k£  are  treated  as  linear  functions  of  temperature 


k  -  k°  +  k*  (T  -  T  ) 
o 


^7 


(9’) 


*ACTO«r  MUTUAL  If  MATCH  CO«K>*ATION 

21011.7 


By  substituting  eq  (8)  in  eq  (1)  and  rearranging  terms,  one  obtains: 

*  3  ^ 


(pc  +  p  c  )  -  ~- 

'^a  pa  Mc  pc'  3t  3x 


k 

s  3x 


+  u  w.V  -  sf  ,{<!  +  (V  hc 


(10) 


This,  rather  than  that  of  eq  (1),  is  the  form  of  the  energy  conservation  equation 
which  is  actually  solved  by  the  model.  The  convective  term  was  not  expanded  in 
order  to  keep  the  finite-difference  scheme  conservative. 

The  two  components  of  the  convective  term  are: 


3h 


M 


g  3x 


and 


3M 
i  — 8 
g  3x 


^g  3t 


The  first  represents  a  volumetric  source  or  sink  of  energy  for  the  solid,  due  to 

the  fact  that  the  specific  enthalpy  of  the  gas  mass  flow  M  is  changing.  The 

8 

second  component  identifies  the  existence  of  a  net  energy  loss  from  the  solid 
equal  to  the  sensible  energy  of  the  gases  produced  in  the  pyrolysis  reaction. 

When  this  latter  component  is  combined  with  the  last  term  in  eq  (10) ,  the  factor 
multiplying  3pg/3t  becomes: 

Q*  -  Q  +  (hfl-  hc*Pf/Pl)/(l  -  Pf/Pi)  -  h  (11) 

Kung,  in  his  paper*  (see  footnote  to  first  page  of  this  Appendix),  Illustrates 

* 

the  relevance  of  this  temperature-dependent  heat  of  pyrolysis,  Q  .  His  deriva¬ 
tion  is  repeated  here  for  two  reasons: 

1)  Values  for  the  heat  of  pyrolysis  quoted  in  the  literature  are  often  values 

* 

for  Q  and  not  for  the  constant  heat  of  pyrolysis  Q  at  reference  temperature 

Tq  (at  T«Tq,  hfl  ■  h£  ■  hg  ■  0) ,  used  here.  As  a  consequence,  proper  care  should 

be  taken  in  adopting  values  for  the  heat  of  pyrolysis  recommended  by  other 

2 

authors  (the  problem  has  also  been  discussed  by  this  writer  elsewhere  ) . 

2)  Users  of  subroutine  SPYVAP  may  be  tempted  to  delete  from  the  model 

the  volatiles  -  solid  heat  exchange  by  setting  c  >0.  They  should  realize, 

P8 
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however,  chat  this  would  have  the  net  effect  of  also  changing  the  energetics  of 

the  pyrolysis  reaction  through  the  disappearance  of  h  from  eq  (9) .  In  the  cur- 

8 

rent  version  of  SPYVAP  it  is  not  possible  to  eliminate  the  first  component  of 
the  convective  term  without,  at  the  same  time,  canceling  the  second. 

The  general  form  of  the  boundary  condition  at  the  front  surface  of  the  slab 


x.0-^,i  +  V(T.,rT.J 


where  q"  net  radiative  flux  received  by  surface, 
h*  convective  heat  transfer  coefficient, 

temperature  of  the  gases  in  front  of  the  surface, 

Tg  surface  temperature, 
c  surface  emissivity, 

8  2  4 

0  Stefan-Boltzmann  constant  (=  5.669*10  W/m  °K  ) 

The  suffix  1  indicates  that  the  quantities  refer  to  the  front  surface.  Since  the 
model  assumes  a  sign  convention  for  the  heat  fluxes,  according  to  which  the  heat 
flux  is  positive  when  entering  the  slab,  the  analogue  of  eq  (12)  for  the  back 
surface  (suffix  2)  requires  a  +  instead  of  a  -  in  front  of  the  conduction  term 
on  the  left  hand  side.  The  quantities  h*  and  e  are  given  constant  values  which 
may  include  0. 

The  net  radiative  flux  received  by  the  surface  is  the  quantity  to  be 
assigned  the  prescribed  time-dependent  values  of  heat  flux,  when  the  problem  re¬ 
quires  a  surface  flux  boundary  condition.  When  the  surface  temperature  Tg  is  pre 
scribed,  the  model  computes  the  heat  flux  implied  by  the  surface  temperature  vari 

ation  and, there fore,  q"  replaces  T  as  an  output  of  the  calculation.  The  possi- 

r  8 

bllity  of  assigning  different  values  to  the  constants  h  and  e  at  both  surfaces 

add6  to  the  flexibility  of  the  model.  As  an  example,  the  case  of  convective 

heating  with  the  front  surface  exposed  to  an  environment  at  temperature  T^  ^ 

is  handled  by  setting  e.  ■  q"  .*•  0  and  by  assign!.  ,  to  h*  the  value  of  the  con- 

1  f)!  1 

vective  heat  transfer  coefficient. 
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A. 3  FINITE  DIFFERENCES  FORM  OF  THE  EQUATIONS 

The  slab  is  divided  in  N  slices  of  constant  thickness  Ax  *  Z/ N.  As  a  result, 

values  of  the  dependent  variables  (T,  p  ,  M  )  are  calculated  at  N+l  equally  spaced 

®  6 

grid  points.  It  is  conventionally  assumed  that  grid  point  1  lies  on  the  front 
surface,  grid  point  N+l  on  the  back  surface.  In  the  interior  of  the  slab  energy 
conservation  is  enforced  for  a  slice  bounded  by  planes  half  way  between  succes¬ 
sive  pairs  of  grid  points .  The  temperature  and  density  are  assumed  to  be  constant 
across  the  slice  and  equal  to  the  values  associated  with  the  grid  point  at  the 
center  of  the  slice.  Near  the  two  surfaces  an  energy  balance  is  Imposed  for  the 
half  slice  extending  from  the  surface  to  Ax/2  below  it.  The  surface  values  of 

T  and  p  are  taken  to  characterize  the  whole  half  slice, 
s 

The  finite  differences  formula  is  obtained  by  integrating  eq  (10)  across  a 
slice  and  by  using  the  Crank -Nicolson  method  to  express  the  different  derivatives. 

At  an  interior  point  i  (1  <  i  <  N  +  1) ,  conservation  of  energy  in  the  step 
from  time  t^  to  «  t^  +  At  is  written  as: 


'  j+Js 

5+k  j+k) 

ks,i+is 

V 

T  -  T 

i+1  Ai 

1+k 
Ca ,  i-J* 


j+4  j+4 

T  -  T 
i-1  i 


1  /  j+ 

)/  S. 


j+%  j+5s  j+*s  j+^s 


i+4  hgti-^  ~  bg.i-^s 


(13) 


-  Ax 


ef  ( 


Q  +  (hfl-  hc«pf/pi)/(l  -  pf/p 


'I  J+*S 


where : 


flfel 

[St  J 


Ps,i  Pf 

~  aP  i’-  p'77.  exp 


-  E/RT 


n 


(14) 
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j+*5  j+*5  r3p 


g»i_Js  g.i-^s 


R) 


J+*s 


•Ax 


(15) 


and  the  mean  enthalpies  are  calculated  from  eqs  (2)  and  (2*). 

Following  a  commonly  used  convention,  subscripts  refer  to  grid  point  location, 
superscripts  to  time  step.  Fractional  subscripts  indicate  location  of  the  cell 
boundary,  i.e.,  mid-point  between  adjacent  grid  points.  Similarly,  fractional 
superscripts  (j+ij)  indicate  average  between  present  (j)  and  new  or  updated 
values  (j+1) . 

With  these  rules  in  mind,  eq  (13)  can  be  written  as: 


j+1 


j+1 


j+1 


(-Vi +  ci-i>  Ti-i +  <Ai -  Bi +  Bi-r  ci  *  c,-i>  Ti  *<-Br  ci>  Vi 


(Bi-i'  ci-i)  Ti-i  +  tAi_  v  Bi-i+  ci~  ct-i>  Ti  + 


(16) 


+  (Bi  4  Cl(  Tl+1  1  4  To  (-C1  4  Cw)  4  0l 


where 


j+^S 

Ai  “  (PaCpa  +  Vpc^ 
J+^5 

Bi  *  kS ,i+%  '  (2,Ax) 


1  *+ 

C  -  7  M  . 
i  A  g,i++i 


o  ,  1  * 

c  _  +  7  C_ 


j+*S 


Dt  -  -  Ax 


pg  2  pf>  [  i+4 

'J*S 


-  T 


/•3p  f  'I 

if  '  h  +  (ha"  hc  pf/pi)/(1  •  pf/pi> 
i 


j+*5 

i 


(17) 

(18) 

(19) 

(20) 


The  form  of  eq  (16)  can  be  further  simplified  by  writing: 


■i  Ti_i  +  ui  Ti  +  T 


j+1 


j+1 


j+1 

i+1  “  bi 


(21) 
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where  the  constants  mit  u^,  and  represent  the  three  coefficients  in  brackets 
in  the  LHS  and  the  whole  term  in  the  RHS  of  eq  (16)  respectively. 

At  a  point  on  the  surface  (i«l  or  i“N+l) ,  the  analogue  of  eq  (13)  can  be 
obtained  by  enforcing  conservation  of  energy  for  the  half  slice  extending  from 
the  surface  to  Ax/2  below  it.  In  the  case  of  the  front  surface  the  energy  balance 
becomes : 


1  ,  .  &  Ax 

T  (P  c  „  +  P  C  ),  7T 


j+i  j 
T.  -  T 


2  a  pa  c  pc  1  At  (  1 

j+b  j4^  j4^  , 

+  M  ,ch  1C-M  .  h  .  -  Ax 
g.1-5  g,1.5  g,l  g.l  2 


3p 


*8,1.5 


j+^5 

r  -  t 

2  *1  J 


/ 


Ax  + 


,3t  j 


Q  +  (h  -  h  p  /p  )/(l-p ,/p  ) 
acil  t  l 


+  q"  ,  +  h  T  -  T.  . 

4r ,1  i  (  “,1  1  J 


i4^)  f  3+^3  j+*s 

£10  (T1  J  ** 1 


(13’) 


where : 


j+*s 

M  .  -  M  .  - 
g.l  g.l. 5 


f^sl 

[3t 


J+^s 


Ax/2 


(13’) 


The  analogue  of  eq  (16)  now  is  : 


(t  V  V  V  V  c„K+1 


-  Er  ci 


1  J  -f 

1 

T.  +  1 

B,+  C, 

J  1  l 

1  lj 

1  3 
T2  + 


(16’) 


+  4  T 


-  C.+  ~  C 
1  2  o 


+  4  D,  +  h*  T  ,  +  q" 

2  1  1  °°,1  r,l 
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where: 


Finally : 


J+l 


*  + 1 

1  +  2 

Cl° 

J+H 

f  o 

'g.l  * 

[Cpg 

j+l 

',‘T 


(18’) 


(19') 


U1  T1  +  fl  T2 


b.  +  4"  n 
1  r  ,1 


(21') 


The  conservation  of  energy  at  the  back  surface  is  written  in  a  similar 
way,  leading  to: 


where: 


Finally: 


j+l 


<-v  V  T»  + 


Aj^+1+  B„  . ,  +  B„-  Cv, , ,  +  C. 


N+l  N  N+l  Nj 


3  J+l 

T. 


N+l 


(B  -  C  )  T  + 
'  N  V  N 


2  VkT  Vf  BN  +  CN+l"  CN 


T  + 
N+l 


+  4  T 


[Vf  CN+lJ  +!Vl+h2T-,2  +  qr,2 


®N+1  “  2  h2  +  2  C2° 


i  ^ 

CN+1  “  2  Mg ,N+1 


f  o  1  *  f, 
[CP8  2  Cpg[ 


j+*s  ) 


T  -  T 
N+l  o 


j+l  J+l 

“n+1  XN  +  “n+I^+I  “  bN+l+  ^r  ,2 
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Equations  (21*),  (21")  and  eq  (21)  for  i«2,N  form  a  system  of  N+l  equations 

in  N+l  unknowns.  The  unknown  quantities  are:  for  i*2,N  and  T^+\ 

(when  surface  boundary  condition  is  on  flux)  or  q"  1 ,  q"_  (when  surface  boundary 

r ,  i.  r9l 

condition  is  on  temperature).  The  solution  of  the  system  of  equations  is  found 

by  using  an  algorithm  for  the  inversion  of  tridiagonal  matrices. 

The  switch  from  a  surface  boundary  condition  where  q^j  is  prescribed  to  a 

surface  temperature  condition  (where  T,+^  or  is  prescribed)  is  accomplished 

^  1+1  1+1 

by  interchanging  the  positions  of  q"  .  (or  q"  _)  and  T^  (or  Tj;  ,)in  the  set 

XT  y  X  r  j  a  X  N't  X 

of  equations.  For  example,  when  the  heat  flux  q^  ^  is  prescribed  at  the  front 
surface,  the  first  two  equations  of  the  system  read: 


J+l 


j+l 


U1T1 


+  f  T 
rl  l2 


V  K,i 


(22) 


j+l 


j+l 


j+l 


m2  ^1 


+  u2  T2 


+  f  T 
2  3 


(23) 


In  the  above  equations,  all  the  unknown  equatities  appear  to  the  left,  all  those 
that  are  known  to  the  right  of  the  equal  sign.  If  the  problem  prescribes  a 
surface  temperature  history,  then  is  known  and  q^  ^  is  the  quantity  to  be 

determined.  The  positioning  of  known  and  unknown  variables  on  separate  sides 

of  the  equal  sign  can  then  be  enforced  by  interchanging  u  T^+^  and  q"  in 

eq  (22)  and  by  replacing  eq  (23)  with  the  sum  of  itself  and  eq  (22) .  Upon 

reordering  of  the  terms  in  the  second  equation,  the  alternate  set  finally  is: 


-  K.i-  fiTi+1 


j+i  j+i 


-  qjfl+  <fi+  u2)  T2  +  f2T3 


j+l 


V  b2  "  (ul+  m2}  T1 


(22’) 


(23') 


In  a  similar  way,  the  last  two  equations  of  the  system  can  be  modified  to  allow 
for  a  surface  temperature  condition  at  the  back  surface. 
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As  mentioned  earlier,  the  model  offers  the  option  of  prescribing  the  boundary 
condition  below  the  surface,  rather  than  at  the  surface.  This  is  done  through  a 
two-step  iterative  procedure,  which  selects  the  surface  temperature  value  that 
satisfies  the  prescribed  "below  the  surface"  temperature  history.  The  procedure 
and  the  provisions  to  delay  instabilities  are  further  Illustrated  in  the  section 
describing  the  details  of  subroutine  SPYVAP.  However,  the  potential  user  should 
be  aware  of  the  fact  that  those  provisions  may  not  be  sufficient  in  situations 
where  the  temperature  condition  is  prescribed  at  a  location  too  far  from  the 
surface.  A  more  precise  definition  of  what  depth  is  too  far  will  require 
preliminary  numerical  experimentation  on  the  part  of  the  user  for  the  particular 
case  to  which  the  model  is  being  applied. 
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A. A  COMPUTER  PROGRAM 
A. 4.1  MAIN  Program 

The  MAIN  program,  whose  listing  is  given  in  one  of  the  following  sections, 
was  designed  for  a  particular  application,  in  which  some  of  the  input  data  were 
entered  through  cards  while  others  were  read  from  a  disk  file.  Clearly,  dif¬ 
ferent  users  will  have  different  requirements  and  they  will  need  to  modify 
MAIN  accordingly.  Despite  its  lack  of  generality,  this  part  of  the  program  is 
reported  here  to  illustrate  what  variables  require  to  be  initialized.  Under¬ 
standing  the  current  version  of  MAIN  and  the  format  of  the  output  (subroutine 
OUTPUT)  is  all  that  is  required  in  order  to  use  the  computer  model. 

A  list  of  the  FORTRAN  variables  used  in  subroutime  SPYVAP  is  given  in 
the  next  section.  For  convenience,  those  which  must  be  initialized  in  MAIN 
are  also  reported  here. 

a)  Grid  Geometry 

DX  slab  thickness,  £. [m] 

N  number  of  slices  (number  of  grid  points  ■  N  +  1) 

b)  Thermal  Properties  and  Pyrolysis 

CPA,  CPC,  CPG  specific  heats  of  active  material,  char  and  volatiles, 

Cpa’  cpc’  Cpg  Uoules/Kg'C]  (eq  2') 

CPAS,  CPCS,  CPGS  temperature  coefficients  of  specific  heats,  c _ , 

DARCY 

PF 

QP0 

RHOF 


pa 


c  ,  c  [joules/Kg*C2]  (eq  2') 
pc  pg 

parameter  for  control  of  migration  of  volatiles 
(■  1.  for  condition  of  eq  (4);  f  1.  for  flow  in  the 
direction  of  decreasing  density) 
pre-exponential  factor,  a 


[s'1]  (eq 


7) 


heat  of  pyrolysis  at  reference  temperature  (T  ) , 

o 

Q  [joules/Kg]  (eq  1) 

fraction  of  Initial  density  at  completion  of  pyrolysis, 

pf/pi  M  (eq  5) 

initial  density,  p  [Kg/m^] 


RHOW 
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c) 


d) 


TCA,  TCC 


TCAS,  TCCS 


TRNEG 

T0K 


thermal  conductivities  of  active  material  and  char, 

k°  k°  [W/m  *C]  (eq  9’) 
o  c 

temperature  coefficients  of  thermal  conductivities, 

k*  k*  [W/m  #C2]  (eq  9') 
a  c 

activation  energy,  E  [ joules /Kg-mole]  (eq  7) 

initial  temperature  (also  reference  temperature),  T  [°K] 


Boundary  Conditions 

(J  -  1/2  or  1/2  in  the  variable  name  indicates  quantity  referring  to  the 
front /back  surface) 

EPS1,  EPS2  surface  emissivities ,  (®9  12) 

HI,  H2  convective  heat  transfer  coefficients, 

h!  [W/m2°C]  (eq  12) 


KBC(J) 

TBC(J,1) 
TIMEBC  (I) 
TINF1,  TINF2 

XBC  (J) 


1’  “2 

index  for  type  of  boundary  condition 

(-  1,  temperature  B.C.;  «  2  heat  flux  B.C.) 

2 

boundary  condition  values,  T  or  q^j  [®K  or  W/m  ] 
times  corresponding  to  B.C.  values,  t [ s ] 
temperature  of  gases  in  front  of  slab  surface, 

r.,i-  i.,2  f«  <««  12) 


distance  below  surface  at  which  temperature  B.C.  is 
prescribed,  [m] 

Control  of  Step  Sire  and  Accuracy 


DTIMAX 

DTIMIN 

DTSTEP 

1STEP 

ITERID 

ITERMX 

LASTEP 

EKRMX 

RELAX 


maximum  time  step  [s] 
minimum  time  step  [s] 

desired  mean  temperature  variation  per  Integration 
step  [ °K) 
step  number 

desired  maximum  number  of  iterations  per  integration  step 
absolute  maximum  number  of  iterations  per  integration  step 
maximum  total  number  of  steps 

maximum  error  between  temperatures  from  successive  iterations 
relaxation  of  temperature  values  from  successive  iterations 
(■  0,  no  relaxation;  *  1,  maximum  effect) 
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SLIMIT  maximum  Increase  of  rate  of  change  of  boundary  condition 

value;  expressed  as  a  fraction  of  the  rate  at  previous  step 
(when  SLIMIT  >  0)  or  assuming  | SLIMIT | *DTSTEP  as  limit  on 
mean  temperature  variation  due  to  correction  (when  SLIMIT  <  0 ); 
operates  when  temperature  beneath  the  surface  is  prescribed  . 

TMLAST  time  limit  for  performance  of  integration  (limit  on  physical 

time  and  not  on  computer  time),  [s] 

e)  Output  Control 

NPROF  number  of  integration  steps  between  outputs  of  profile 

variables 

NSTAT  number  of  integration  steps  between  outputs  of  station 

variables 

As  a  general  rule  variables  are  entered  in  their  dimensional  form,  with 

dimensions  in  S.I.  units  (Kg,  m,  s  and  derived  units).  Temperatures  are  in 

degrees  Kelvin  (°K) .  The  variable  1STEP  must  be  initialized  to  0  by  MAIN. 

Note  that  the  pyrolysis  option  can  be  deleted  from  the  model  by  assigning  0  to 

the  pre-exponential  factor  a^  (PF) .  In  that  case  the  thermal  properties  of 

the  solid  are  equal  to  those  of  the  active  material  (c  ,  k  )  at  all  times. 
n  pa  a 

In  the  version  of  MAIN  presented  here  a  large  part  of  the  program  is 
occupied  by  a  set  of  instructions  which  perform  the  operation  of  reading  from 
a  disk  file  the  boundary  values  for  temperature  or  heat  flux  at  the  two  surfaces 
and  the  times  corresponding  to  these  values.  As  a  result,  arrays  TBC(J,I)  and 
TIMEBC(I)  are  initialized;  the  operation  is  controlled  by  the  values  of  the 
parameters  CHF,  CHB,  DBCMN  and  TMX.  The  reader  does  not  need  to  worry  about 
this  section  since  he/she  will  have  to  rewrite  it  in  any  case. 

The  second  function  performed  by  MAIN  is  to  print  the  initial  conditions 
and  a  table  of  the  boundary  values.  The  only  reminder  here  is  that  the  dif¬ 
ferent  quantities  are  printed  in  S.I.  units.  Finally,  a  small  section  is 
dedicated  to  the  supervision  of  the  calls  of  subroutines  SPYVAP  and  OUTPUT. 

The  number  of  Integration  steps  to  be  performed  before  returning  to  MAIN  (IOUT) 
is  the  only  variable  required  by  SPYVAP  in  the  arguments  list.  Similarly,  sub¬ 
routine  OUTPUT  requires  the  argument  IPRINT  to  be  initialized  to  1  for  output 
of  station  variables,  to  2  for  output  of  profiles.  When  IPRINT  -  0  no  output 
takes  place.  Note  that  the  value  of  IPRINT  is  determined  using  the  auxiliary 
real  variables  ANSTAT  and  ANPROF. 
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Execution  terminates  when  any  of  the  following  three  conditions  is  verified: 

1)  current  time  equals  TMLAST, 

2)  total  number  of  integration  steps  equals  LASTEP,  or 

3)  termination  code  IF1N  is  returned  from  SPYVAP  at  a  value  f  0. 

Which  one  of  these  conditions  terminated  execution  of  the  program  can  be 
determined  from  the  last  line  of  the  printed  output. 

A. 4. 2  Subroutine  SPYVAP 

The  only  variable  passed  to  this  subroutine  through  the  argument  list  is 
IOUT,  which  is  the  number  of  steps  to  be  performed  before  returning  control  to 
MAIN.  All  the  other  transfers  are  implicitly  accomplished  by  using  the  common 
area  SPYCOM  for  storage.  Subroutine  SPYVAP  is  organized  in  8  chapters,  the 
first  3  of  which  are  executed  only  when  IST£P«0,  usually  corresponding  to  the 
first  call  of  the  subroutine.  The  different  chapters  are  now  described  in  detail. 
Chapter  0 

The  number  of  the  last  step  to  be  completed  before  returning  to  MAIN  (ISTEPR) 
is  evaluated.  Control  is  then  transferred  to  the  beginning  of  chapter  1  if 
ISTEP  is  equal  to  0,  otherwise  execution  proceeds  from  the  beginning  of  chapter  3. 
Chapter  1 

Current  time  (TIME)  is  initialized  with  the  time  at  which  the  first  boundary 
values  are  available  and  the  time  step  is  set  equal  to  DTIMIN.  The  thickness 
of  the  slice,  Ax,  is  substituted  for  the  slab  thickness,  £,  in  DX.  The  variable 
IBC(J),  with  J«1  for  front  and  J»2  for  back  surface,  represents  the  grid  point 
immediately  to  the  left  of  the  location  where  a  temperature  boundary  condition 
is  prescribed.  IBC(J)~0  when  the  condition  is  at  the  surface.  The  distance 
between  the  interior  boundary  point  and  the  grid  point  immediately  to  its  left, 
expressed  as  a  fraction  of  Ax,  is  stored  in  XBC(J).  A  minimum  temperature  (TL) 
is  determined  below  which  the  pyrolysis  calculation  is  not  performed:  such 
temperature  is  the  temperature  at  which  1  percent  of  the  pyrolyzable  material 
would  be  gasified  for  a  pyrolysis  duration  equal  to  the  heat  diffusion  time 
associated  with  the  slab. 
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Finally  all  specific  heats  are  normalized  with  respect  to  c  ,  which  is 

opa 

stored  in  CPW,  and  all  thermal  conductivities  with  respect  to  k  ,  which  is 
stored  in  TCW. 

Chapter  2 

Temperatures  in  arrays  T(I)  and  TP (I)  are  set  equal  to  the  initial  tempera¬ 
ture  (T0K) .  The  value  1  is  put  in  array  RHO(I),  which  contains  local  den¬ 
sities  normalized  with  respect  to  the  initial  density,  p^/p^.  Other  arrays  and 
auxiliary  variables  are  also  initialized  in  this  chapter. 

Another  quantity  defined  here  is  the  initial  time  rate  of  increase  of  the 
variable  to  be  used  for  the  surface  boundary  condition:  this  quantity  is 
stored  in  array  SLOPE  (J)  (J»l,  front;  J-2,  back  surface).  SLOPE(J)  is  used 
in  the  case  in  which  a  temperature-time  history  at  distance  below  the  sur¬ 
face  is  specified  for  the  boundary  condition.  As  indicated  by  Carlslaw  and 
Jaeger  ("Conduction  of  Heat  in  Solids",  Oxford  University  Press,  p.  388,  1959), 
the  temperature  increase  T  -  Tq  at  a  distance  x  from  the  surface  caused  by  a 
linear  increase  of  surface  temperature  from  Tq  to  Tg  in  time  At  is  given  by: 


Tq  -  (Tfi-  Tq)  (l+2y  )  erf c  (y) 


where : 


2/a*At 


(24) 


(25) 


and  a  is  the  thermal  diffusivity  of  the  material.  The  above  fact  is  recognized 
in  the  program  and  eq  (24)  is  used  by  determining  a  from  the  thermal  properties 
of  the  active  material  at  reference  temperature  and  by  setting  At  equal  to 
DTIMAX. 

Chapter  3 

This  chapter  is  concerned  with  determining  the  boundary  values  to  be  used 
in  the  integration  step  to  be  performed.  This  is  done  by  interpolating  linearly 
the  values  contained  in  array  T2C(J,I).  Note  that  temperature  boundary  values 
are  calculated  at  the  new  time  t^+1  («  t^  +  At) ,  while  heat  flux  boundary 
values  are  evaluated  at  the  mean  between  current  and  new  time  t^**5  (■  t^+  ^  At) . 
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When  a  temperature  condition  Is  prescribed  not  at  the  surface  but  below  It , 

It  Is  necessary  to  determine  the  appropriate  value  for  the  surface  temperature. 

Such  value  Is  calculated  through  a  two-step  Iterative  procedure.  A  first  ap¬ 
proximation  to  the  value  of  the  surface  temperature  is  found  by  assuming  a 
time  rate  of  Increase  equal  to  that  of  the  previous  step.  Then  the  temperature 
profile  in  the  slab  is  calculated  by  solving  the  energy  equation  and  the  value 
at  the  depth,  where  the  temperature  is  prescribed,  is  compared  with  the  boundary 
value . 

On  the  basis  of  this  comparison  a  second  approximation  to  the  surface  tem¬ 
perature  is  found.  In  doing  so,  the  scaling  factor  implied  by  eq  (24)  is 
properly  taken  into  account,  by  introducing  approximate  expressions  for  the 
error  function.  A  new  temperature  profile  is  found  and  a  third,  final  value 
for  the  surface  temperature  is  calculated  by  linear  interpolation  between  the 
first  two  surface  temperatures  and  the  corresponding  temperatures  at  the  pre¬ 
scribed  depth. 

In  general,  the  procedure  will  work  smoothly  for  a6  long  as  the  value  of 
the  quantity  y  (eq  (25))  is  not  too  large.  A  limit  on  the  per-step  variation 
of  the  rate  of  change  of  the  surface  temperature  is  imposed  through  the  variable 
SLIMIT  (SLIMIT  ■  1  allows  a  maximum  variation  of  100  percent  on  SLOPE  (J)  when 
y  ■  0) .  However,  in  some  cases  a  considerable  amount  of  numerical  experimenta¬ 
tion  with  different  values  of  the  control  parameters  may  be  necessary  to  avoid 
the  onset  of  oscillations  in  the  values  of  the  surface  temperature. 

Chapter  4 

After  the  appropriate  set  of  boundary  conditions  has  been  found,  this 
chapter  becomes  the  top  of  the  iteration  loop.  First,  density  changes, 

DRHO(I),  are  calculated  from  eq  (14)  and  the  mean  densities  during  the  time 
step  p  J  /p  are  stored  in  array  RHOA(I) .  Then,  the  updated  distribution  of 

S  ,  X  X 

the  flux  of  volatiles  is  calculated  from  eq  (15):  when  DARCY  ■  1.  the  con¬ 
dition  of  eq  (4)  is  implemented;  otherwise  the  volatiles  are  assumed  to  flow 
in  the  direction  of  decreasing  densities.  Mean  values  of  thermal  conductivity, 

k  and  specific  heat,  are  calculated  from  eqs  (9)  and  (8’)  respectively. 

ps,i 
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The  convective  term  and  the  energy  source  due  to  pyrolysis  are  then  evaluated. 

By  the  end  of  this  chapter,  the  quantities  A^,  B^,  and  defined  in 
eqs  (17) -(20)  have  all  been  assigned  their  respective  values.  Because  of  the 
order  in  which  the  different  operations  are  performed  by  the  program,  the  in¬ 
structions,  whose  execution  is  superfluous  wher.  some  of  the  options  of  the  pro¬ 
cedure  are  not  in  use,  are  simply  bypassed  with  a  consequent  saving  in  execution 
time. 

Chapter  5 

All  the  operations  relating  to  the  tridiagonal  matrix  are  performed  here. 

the  coefficients  m^,  u^,  f  and  b^  in  eqs  (21),  (21')  and  (21")  are  evaluated 

first;  when  the  boundary  condition  is  on  the  temperature  of  the  surface,  the 

values  of  some  of  the  coefficients  are  modified  as  implied,  for  example,  by 

eqs  (22')  and  (23')  in  the  case  of  the  front  surface.  The  tridiagonal  matrix 

is  then  inverted  using  the  algorithm  discussed  by  Forsythe  and  Moler  ("Computer 

Solution  of  Linear  Algebraic  System,"  Prentice-Hall,  p.  118)  and  the  solution 

values  are  stored  in  BT(I).  Taking  again  the  front  surface  as  an  example, 

it  is  realized  that  when  the  boundary  condition  is  on  the  surface  temperature, 

BT(1)  contains  q"  and  not  T^+^  and,  therefore,  the  appropriate  switch  is 
r»-L 

made  so  that  BT(I)  contains  the  temperatures  T^ 

Chapter  6 

This  part  of  the  subroutine  performs  the  following  operations: 

1)  update  of  TP(X)  array  containing  the  temperatures  at  the  new 

time  and  check  of  the  difference  with  respect  to  the  result  of  the 
previous  iteration,  with  decision  on  whether  to  perform  another 
iteration  or  accept  the  calculated  profile  as  sufficiently  accurate; 

2)  evaluation  of  new  step  size  to  be  used  in  the  next  step  and  deter¬ 
mination  of  extrapolated  values  for  temperatures  T^+^  and  dimension¬ 
less  density  increments  for  the  next  step; 

3)  evaluation  of  total  energy  that  has  entered  the  slab  through  the 
front  (Q1DT)  and  the  back  surface  (Q2DT) . 
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A  few  additional  comments  will  Illustrate  the  operations  carried  out  In 
the  chapter.  When  the  accuracy  requirement  is  not  met,  the  option  is  available 
to  relax  the  temperature  profile  using  the  profile  from  the  previous  iteration. 
The  maximum  weight  assigned  to  the  temperatures  from  the  previous  iteration  is 
equal  to  the  value  of  the  variable  RELAX.  The  value  of  such  weight  is  decreased 
automatically  as  the  residual  error  approaches  the  maximum  allowed  (ERRMX) . 

A  series  of  controls  are  available  to  optimize  the  step  size.  First  of  all, 
whenever  the  number  of  required  iterations  reaches  ITERMX,  the  step  size  is 
halved,  a  message  is  printed  and  the  calculation  is  restarted  from  the  current 
time.  Upon  successful  completion  of  a  step,  the  value  of  the  step  size 
(DTIME)  for  the  following  integration  is  evaluated  by  insuring  that  the  average 
of  the  absolute  temperature  variation  in  the  slab  remains  close  to  DTSTEP. 

DTIME  is  decreased  if  the  integration  Just  completed  required  a  number  of  iter¬ 
ations  greater  than  ITERID.  Finally,  execution  stops  when  the  value  calculated 
for  DTIME  is  lower  than  DTIMIN. 

Chapter  7 

The  step  number  counter  (ISTEP)  is  updated  and  a  decision  is  made  whether 
to  return  control  to  MAIN  or  go  to  the  beginning  of  chapter  3  for  another  in¬ 
tegration.  Return  to  MAIN  can  be  caused  by  any  of  the  following  conditions: 

1)  current  time  is  greater  than  or  equal  to  IMLAST, 

2)  step  counter  shows  that  a  number  of  steps  equal  to  LASTEP  has  been 
completed , 

3)  a  return  code  other  than  0  is  in  IFIN,  or 

4)  the  number  of  steps  (I0UT)  required  for  the  current  call  of  the  sub¬ 
routine  has  been  completed. 

A. 4. 3  Subroutine  OUTPUT 

This  subroutine  is  dedicated  to  the  output  of  station  variables  and  pro¬ 
files  with  the  frequency  implied  by  the  values  of  NSTAT  and  NPROF  in  MAIN.  Pro¬ 
files  are  printed  only  when  the  parameter  IPRINT  is  greater  than  1;  in  addition, 
the  density  and  mass  flux  values  (at  the  grid  points)  are  printed  only  when 
pyrolysis  is  taking  place  as  evidenced  by  a  value  of  RHOBAR  less  than  1. 
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Note  that  the  quantities  appearing  In  the  output  are  all  associated  with  the 
nean  time  t^+  y  At  (TIMEAV)  :  the  arrays  printed  for  the  temperature 

and  density  profiles  are  TBAR(I)  and  RHOA(I)  respectively. 

The  only  exception  is  in  the  additional  output  referring  to  the  boundary 
temperature  match  :  this  is  active  when  a  temperature  boundary  condition  is 
prescribed  below  the  surface  instead  of  at  the  surface.  The  temperatures  shown 
are  the  value  that  the  boundary  temperature  should  have  at  time  t^+^*  t^+  At 
and  that  implied  at  the  same  time  by  the  computed  profile  (T(I)).  The  differ¬ 
ence  between  the  two  values  gives  an  indication  of  how  closely  the  temperature 
boundary  condition  is  being  enforced  by  the  model. 

The  three  components  (see  eq  (12))  of  the  net  surface  heat  flux  are  re¬ 
ported  among  the  station  variables.  They  are: 

1)  radiative  flux  q£,  (QRAD1,  QRAD2) , 

2)  convective  flux  h*^-  T^)  ,  (QC0NV1,  QC0NV2)  , 

3)  surface  reradiation  eoT^  ,  (RERAD1,  RERAD2) . 

s 

The  integrals  for  the  two  surfaces  of  the  net  fluxes  are  given  in  Q1DT  and  Q2DT. 
The  quantities  MG1  and  MG2  are  the  blowing  rates  at  the  two  surfaces,  positive 
when  out.  As  a  reminder,  all  dimensional  quantities  are  shown  in  S.I.  units. 
More  specifically: 

times  in  s 

mass  fluxes  in  Kg/m^s 

2 

energy  fluxes  in  W/m 

temperatures  in  °K. 
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A. 5  LISTING  OF  VARIABLES  IN  SUBROUTINE  SPYVAP 


The  meaning  of  the  FORTRAN  variables  which  appear  in  the  subroutine  SPYVAP 
is  given  in  the  following  list.  In  order  to  facilitate  the  understanding  of 
the  subroutine  the  explanation  of  the  different  variables  is  accompanied  by  a 
cross  reference  listing,  which  indicates  the  numbers  of  the  statements  where 
the  individual  variables  are  used.  As  a  general  rule,  the  index  I  is  used  to 
indicate  grid  points  (1,  NP1  range),  while  J  distinguishes  between  front  (J*l) 
and  back  surface  (J*2) . 

A(I)  coefficient  defined  in  eq  (17) 

B(I)  coefficient  B^  defined  in  eq  (18) 

BC(J)  value  of  boundary  condition  (on  temperature  or  heat  flux) 

BCTRl(J),  BCTR2(J)  temperature  calculated  from  the  first  and  second 

iteration  when  temperature  is  prescribed  at  a  given  distance 
below  the  surface 

BETA  actual  amount  of  relaxation  of  temperature  profiles 

BETAM1  1.  -  BETA 

BT(I)  coefficient  b^  appearing  in  eqs  (21),  (21')  and  (21") 

C(I)  coefficient  defined  in  eq  (19) 

CHG  maximum  absolute  temperature  change  per  step 

CPA,  CPC,  CP G  specific  heats  at  reference  temperature  for  active  material, 

char,  volatile  products,  c°  ,  c°  ,  c° 

*  pa'  pc’  pg 

CPW  ■  CPA  referen*  \  specific  heat 

*  * 

CPAS,  CPCS,  CPGS  temperature  coefficients  of  specific  heats,  cpa9  cpc, 
cpg  (see  eq  (2’)) 

CPGDA  -  .25*CPG 


CPGDA 


DARCY 


ZJ2  (see  eq  (19’)) 

CM/2  (see  eq  (19")) 
coefficient  D^  defined  in  eq  (20) 

parameter  for  type  of  volatiles  migration:  when  *1  condition 
of  eq  (A)  is  satisfied,  otherwise  migration  follows  the  direction 
of  decreasing  densities  of  the  solid 
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DRHO(I) 

DQ1DT,  DQ2DT 

DTIMAX,  DTIMIN 
DTIME 
DTSTEP 
DX 

EPS1,  EPS2 

ERR 

ERRMX 

F(I) 

FACT 

HFLUX1.  HFLUX2 
HI,  H2 

I 

1BC(J) 

ICPL 

IER 

IFIN 

IOUT 

ISTEP 

ISTEPR 

ISURF(J) 


dimensionless  density  change,  P^)/p^ 

increments  of  energy  stored  in  the  slab  through  front  and 

back  surface 

maximum  and  minimum  time  step  size 
size  of  time  step 

desired  mean  temperature  variation  per  step 
spacing  between  grid  points  Ax  (Initially  slab  thickness) . 
product  of  Stefan-Boltzmann  constant  with  emisslvity  at  front 
and  back  surface  (initially  emisslvity) 

error  between  temperature  profiles  from  successive  iterations 
maximum  accepted  value  for  ERR/CHG  or  ERR  (in  ®K) ,  whichever  is 
smaller 

coefficient  f^,  appearing  in  eqs  (21),  (21*)  and  (21") 
scale  factor  used  for  determining  temperature  values  at  the 
surface  when  temperature  is  prescribed  inside  the  slab 
radiative  heat  transfer,  q^,  at  front  and  back  surface  (positive 
when  going  in) 

* 

convective  heat  transfer  coefficient,  h  ,  for  front  and  back 
surface 

grid  point  index 

first  grid  point  to  the  left  of  location  where  temperature 
boundary  condition  is  imposed  (■  0  for  surface  boundary  condition) 
■  N4-2  -  I 

grid  point  of  maximum  temperature  error  between  successive 
iterations 

termination  label  (■  0  for  normal  run) 

number  of  steps  to  be  completed  by  subroutine  before  returning 
to  MAIN 

current  step  number 

last  step  to  be  completed  by  subroutine  before  returning  to  MAIN 
“  1  for  J*l,  ■  N+l  for  J»2 
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ITER 

ITERID 

ITERMX 

J 

K 

KBC(J) 

LASTEP 

M(I) 

MG(I) 

MG1,  MG2 

N 

NP1 

PF 

QP 

QP0 

Q1DT,  Q2DT 

RCDX 

RJDEN 

RELAX 

RHO(I) 

RHOA(I) 

RHOF 

RH0FM1 

RHOW 


number  of  completed  Iterations 

desired  number  of  Iterations  per  step 

maximum  number  of  iterations  per  step 

index  indicating  front  (J*l)  or  back  surface  (J*2) 

index  referring  to  current  boundary  values  in  TBC  (J,K)  and 

TIMEBC(K)  arrays 

■  1,  temperature  boundary  condition; 

■  2,  heat  flux  boundary  condition 
maximum  number  of  integration  steps 

coefficient  m^,  appearing  in  eqs  (21),  (21')  and  (21") 

mass  flux  of  volatiles,  M  (see  eq  (15)) 

8 

mass  fluxes  of  volatiles  at  front  and  back  surface  (positive 

when  out  of  the  slab,  cf  eq  (15') 

number  of  slices  in  slab 

number  of  grid  points  (*  N+l) 

pre-exponential  factor,  a^ 

last  term  in  brackets  in  eq  (10) 

heat  of  pyrolysis  at  reference  temperature,  Q 

energy  stored  in  slab  through  front  and  back  surface 

-  RHOW*CPW*DX 

factor  used  in  determination  of  coefficient  for  relaxation  of 
temperature  profiles 

maximum  value  of  coefficient  for  relaxation  of  temperature 
profiles 

dimensionless  density  at  current  step,  p^/p. 

1+1  1  *  1 

dimensionless  mean  density  (p^  +  P^)/2p^ 
dimensionless  final  density  Pf/P^ 

-  1.  -  RHOF 
initial  density, 
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SL2MIT  maximum  variation  of  the  rate  of  change  of  boundary  values 

when  boundary  condition  is  prescribed  on  temperature  inside 
the  slab 

SLOPE (J)  rate  of  change  of  boundary  values  at  current  time  step 

T(I)  temperature  at  current  time, 

TBAR(I)  mean  temperature,  (T^+^+  T^)/2 

TBC(J,K)  boundary  values  for  temperature  or  heat  flux 

TCA,  TCC  thermal  conductivities  of  active  material  and  char  at  the 

reference  temperature,  k°,  k° 

a  C 

TCAS,  TCCS  temperature  coefficients  of  thermal  conductivities, 

k  ,  k  (see  eq  (9')) 
a  c 

TCW  reference  thermal  conductivity 

TCWDHX  -  . 5*TCW/DX 

THDIFF  -  TCW/RHOW/CPW  reference  thermal  diffusivity 

TIME  current  time,  t^ 

TIMEAV  mean  time,  t^4***  t^+  At/2 

TIMEBC(K)  times  associate  with  boundary  values  in  TBC(J,K) 

TIMEDT  new  time,  tJ+1-  tJ+  At 

TINF1,  TINF2  temperature  of  the  environments  facing  the  front  and  back 
surface  of  the  slab 

TL  minimum  pyrolysis  temperature 

TP (I)  new  temperatures,  i+1 

TRNEG  initially  activation  energy  E,  then  set  to  -  E/R 

TSLCR(J)  new  rate  of  change  of  boundary  values 

TSTRl(J),  TSTR2(J)  surface  temperature  values  from  first  and  second  iteration 
1KLA5T  maximum  time 

T0K  reference  and  initial  temperature 

U(I)  coefficient  u^,  appearing  in  eqs  (21),  (21*)  and  (21") 

XBC(J)  initially  distance  from  the  surface  of  the  boundary  location 

at  which  temperature  is  prescribed,  then  distance  of  same  location 
from  the  first  grid  point  to  the  left  expressed  as  fraction  of  Ax 
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dimensionless  distance,  y  (see  eq  (25)) 

Initial  values  of  XBC(J) 

attenuation  coefficient  for  in-depth  temperature  change 
(see  eq  (24)) 
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FORTRAN  LISTING 
MAIN  Ingram 


SN  0002  C0MM0N/SPVC0M/BC<2»  , CPA, CPC  ,C  PG, CPAS  ,CPCS  ,CPGS  ,  DARCY,  f'T  ]Mtt, 

1  Or|ME,DTrNIN,DTSTEP,DX,FPSl,EPS2,ERRNX,HFlUXitHFlJX2, 

2  H1,H?,IBCI2I ,IFIN,I STEP,ITER,ITER!0,ITERMX,KBC(2), 

3  LASTFP,MG(441 ,MG1 ,MG2,N,NP1 ,PF ,QP0 ,Q10V , 02 0T , REl  AX, 

4  RH0(  441  •  PH0AI  44  I  ,RH0F , RHOM, SL I  Ml T  ,T  (  44  )  ,  T  BAR  ( 44 1  , 

5  TBCI  2,1501  ,TCA,  TCC  ,  TCA  S.TCCS  ,TI  ME  ,T  1  ME  AY  ,T  1  MEBC  US'* )  , 

6  TINE 1,  TT  NF  2  »  TMLAST ,TRNFG»T0K,XBC(2I 


SN 

0003 

PEAL  MG, MG 1 , MG 2 

SN 

0004 

INTEGER  DIAV1 ,DIAV2,0I AV3 

SN 

0005 

DEFINE  FILP  5(  13000,  80, t.IAVl 

r . 

c- 

— 

-  THERMAL  PROPERTIES 

SN 

0006 

110 

CONTINUE 

SN 

0007 

READ! 5, 100,EN0  =  3201  C PA ,CPC ,C PS ,C PAS ,CPCS  ,C  PGS 

SN 

0008 

READ! 5, 1001  TC A  ,TCC , TC A S , TCC S , PHOW 

c- 

— 

-  PYROLYSIS  CONSTANTS  AND  SLAB  THICKNESS 

SN 

0009 

«EAO( 5, 1CC1  OX, RHOF, PE ,TRN€G , OPO ,TOK , DARCY 

c- 

- PAR  A**FTFR  S  FOR  BOUNDARY  CONDITIONS  AND  GRID 

SN 

0C10 

READ! 5, 1051  KBC ( 1 1 ,KBC (2 1 ,N,l A  STEP, NSTAT .NPROF , IT FPMX, I  TER  ID 

SN 

0011 

ISTEP=0 

SN 

C012 

REA 0(5,  1001  Hl,H2,TINFl,T!NF7,tPSl,EPS2,XBCm,XBr(  2) 

SN 

0013 

RFAO( 5, 1001  OTIMAX.OTI MI N.DTSTE P .TMLAST ,ERR MX  ,  PEL  AX , SL I M IT 

SN 

0014 

100 

FORMAT! 8F10.3) 

SN 

0015 

105 

FORMAT! 81 101 

SN 

00]6 

ANSTAT=NSTAT 

SN 

0017 

anppof=nprof 

r  - 

SN 

0018 

R  FAD( 5, 100 )  CHF ,CH8,0BCMN,TMX 

SN 

0019 

CH=CHF 

SN 

00  20 

CH2=CHR 

SN 

0021 

IF(CHF.LT.CH5 1  GO  TO  115 

SN 

00?3 

CH1=CHB 

SN 

0024 

C  H2=CHF 

SN 

0O?5 

115 

CONTINUE 

SN 

0026 

J*3 

SN 

0027 

IE(CHF.EO. 1000.  1  J*1 

SN 

0079 

iF(CHP.eo. looo.i  j*2 

SN 

0031 

I AV* 168 l 

SN 

0032 

OlAVl*lNT(  (CHI-1.  1  t\  0.  1 

SN 

0033 

I  1*  !NT(CH1»-10*D!AV1 

SN 

0034 

0IAV2*-l 

SN 

0035 

I  2*  1 

SN 

0036 

IE( CH2.F0. 1000. )  GO  TO  120 

SN 

0038 

0IAV2=INT(  (CH2-1.)  /10.  1 

SN 

0039 

12=  INT( CH21- 1040IAV2 

SN 

0040 

0  IA  V2*D I A V2-0I A  Vl-1 

SN 

0041 

120 

CONTINUE 

SN 

0042 

0  IA V3*18-DIAV1-0IAV2 

SN 

0043 

1*1 

SN 

0044 

125 

CONTINUE 

SN 

0045 

READ! 8* IAV, 1301  Ml N, SEC 

SN 

0046 

130 

FORMAT! 6X,I 2, 2X,F 5.2,65X1 

90 
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SI  IK»*T 
SM  n«A8 
ISM  00*9 
ISM  0050 
SN  0051 
SM  005? 
SM  005* 
ISM  0055 
f  SN  0056 
SM  0057 
SM  P"58 
SM  0059 
ISM  0060 
SM  006? 
SM  0063 
SM  0"6* 
SM  "065 
SM  0067 
SM  0068 
SN  0070 
SM  0071 
SM  0073 
SM  0076 
SM  0077 
SM  "079 

SM  "OR" 

SM  008? 
SN  0083 
SM  0" 66 
SM  0085 
SM  0086 
SN  0088 
SM  "090 
SN  0091 
SM  0092 
SM  "09* 
SN  0095 


SM  "096 
SM  0097 


SM  0098 
SN  0099 


SM  "100 
SM  0101 

SN  010? 
SM  0103 


T(IKK(II-te.tfUMT(lf«|»irC 

uv>u««8tm 

READ!  8*  IAV, 1351  I Tf IRI ,1 R-l ,1 01 
135  FORMAT! 10F8. 3) 

TBC1 1,1  )«T! li) 

1FID1AV2.E0.-1I  GO  TO  1*0 
IAV*IAV«DIAV2 

READt  8*  IAV,  1351  I  T{  IR )  ,1  R*1  ,101 
1*0  CONTINUE 

IAV=IAV*0IAV3 
F  INDt  8* IAV) 

TBCt 2,1 )*T! 12) 

ICtCHF.LT.CHB)  GO  TO  1*5 
TBCt  2,1  )  =  T3Ct 1,1) 

TBCt 1,1  )=T( 121 
1*5  CONTINUF 

I  Ft  J.EQ. 3)  GO  TO  150 
TBCt J, I  1  =  0. 

I  Ft  KBC I J).E0.1)  TBCt J,I)*T0K-273. 

150  CONTINUE 

IFt  I.FO.l)  GO  TO  155 
IFt TtMEBCf I ) .GT.TMLASTI  GO  TO  160 
IFt TIMEBCt I I.GE .TIMFBC t!-I)*TMX)  GO  TO  165 
I F  f  ABSt  TBCt  1,  D-TBCI 1  , 1  - 1  )  )  .LT.DBCMN)  GO  TO  125 
155  CONTINUE 

IFt I.E0.150)  GO  TO  160 
1*1*1 
GO  TO  125 
16"  CONTINUE 
IR*  I 

IFt  IR.FQ. 15O.AN0.TMLAST.GT. T!  MF8C  t) 53) )  TML AST  =  T I  MB BC fl 50 ) 
IFtKBCt l).E0.2.AM0.KBCt2).F0.2)  GO  TO  170 
DO  165  J*  1, 2 
00  165  1*1, IR 

IFtKBCt Jl.EO.l)  TBCt J,!!=TBCf J,I)*273. 

165  CONTINUE 
170  CONTINUE 
C 

C - OUTPUT  OF  INITIAL  CONDITIONS  - 

WP|TFt6,21P)  N,DX»TOK,PHOW,PHOF,PF  ,TRNEG  »QP3  ,  DARCY 
210  FORMAT!  'lGEOMETRY,  INITIAL  CONDITIONS  AND  PY POLYS  IS:*/*  N=’, 

1  I3,2X,’ST*’,1PE10.3*2X»,TOK*’,F10.3»2X»,RHOW*,,E10.3*2X, 

2  ’RHOF*’ ,E10.3,2X,,PF  *' ,EI 0. 3 , 2 X ,  ’TRNEG*’ , Eli . 3, ?X , • OPO* • , 

3  E10.3/'  DARCY*’ ,0PF3.0) 

WRITE  I  6,220)  CPA,CPAS,CPC,CPCS,CPG,CPGS,TCA,TCAS,TCC,TCCS 
220  FORMAT! ’OThFRMAL  PROPF P TI E S : ’ / •  SPFCIFIC  HF ATS : ’ , 1 8X , • C PA* • , 

1  1 PF10.3,2X,’CPAS*’ , 

2  F10.3.2X, ’CPC*’,F10. 3,2X,’CRCS=  •  ,E10. 3 ,2X , • CPG*’ ,F10.3,?X, 

3  ’CPGS*  *  ,E  10.  3/’  THERMAL  CONDUCTI V! TI ES : • ,1 3X , • TC A*' , E 13 . 3, 2X, 

*  •TCAS*’,F10.3,2X,’TCC*’,E10.3,?X»,TCCS*’ , El  0.3) 

WRITE! 6,230)  KBCtll ,XBCtll , TI NF 1 ,H1 ,FPS 1 ,KBC 1 2 ) , X BC t 2 ) , T INF 2, 

1  H2.EPS2 

230  FORMAT! ’OBOUNDARY  COMO! TI ONS: ’ /2? X , ’KBC • ,8X , • XBC’ , 1 1 X, • T IMF • , 

1  12X, ’H* ,13X, ’EPS’/’  FRONT  SURF  ACE ; • , I  10 , 1 P*E 1 6 . 3/ •  BACK  SURFA 

?’  , Ill* *E 15. 31 

WRITE!  6,  250)  t  TIMEBCU  )  ,1*1  ,1  B) 

250  FORMAT! *0TIME',*X,1Pl)E11.3/t9X,llE11.3)) 
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SN  9104  MBfTE(4,*«0t  CHF.< T0C« 1 ,1 » ,1-1 , l«l 

SN  9109  *00  FORMAT! *OCM«  • •04.0*1011011 .*/•  f«0*T  0C* ,11*11 .9/(91. IIEU .91 » 

SN  0106  MR ITEI 6, 270 )  CHB,ITBCI2,1 » ,I«1 ,1R) 

SN  0107  270  FORMAT! *0CH#  *,F4.0,1P11E11.3/*  BACK  BC  • ,11E11 .3/ !9X, 1 1*11 .3 » » 

SN  0108  MRITEI 6,240)  DTIMAX, DTI  MI N, 0T STEP .ERRMX ,SLI MIT , RELAX,  ITERMX,  HER 

SN  0109  240  FORMAT! *0STEP  CONTROL  PARAMETERS:  •/ •  OT I  MAX** , l PE9. 2. 2X , 

1  *DTTMIN**,E9.2,2X,*0TSTEP**,E9.2,2X,* ERRMX** ,E9.2,2X, 

3  *SLIMIT**,E9. 2,2X, 

2  *RELAX**,E9.2,2X,*ITERNX«*,0PI3,2X,*!TERI0-*,I3/*  1* ) 

C 

C -  CALL  OF  SUBROUTINE  SPYVAP  - 

SN  0110  310  CONTINUE 

SN  0111  IOUT=!ISTEP/NSTATM)*NSTAT 

SN  0112  IOUTsMINO! I0UT.I I STEP/NPROFyI  ) 4NPP0F) 

SN  0113  IOUT*IOUT-I STEP 

SN  0114  CALL  SP WAP !  1  OUT ) 

C - CHECK  FOR  OUTPUT  - 

SN  0115  IPR INT*0 

SN  0116  IFCELOAT! I  STEP /NSTATI . E 0. FLOAT! I STFP) / ANSTAT)  IPR1NT  =  1 

SN  0118  IFIFLOAT!  1 STFP/NPPOF ) . FQ. FLOAT! I  STEP) / ANPROF)  IPRINT*2 

SN  0120  IF!  IFIN.NF.0.0P.T1MF.GE.TMLAST.0R.ISTEP.F0.LASTEP)  1  PR  I NT=  2 

SN  0122  IF!  IPRINT.GT.Ot  CALL  OUTPUT! I PP I  NT ) 

SN  0124  IF! IFIN.EQ.C.ANO.TIME.LT. TML AST . AND. ISTEP.NE. LAST EP)  GO  To  810 

SN  0126  MPITE16.200)  l STF P ,LASTEP , T1 ME , TML AST , I F 1 N 

SN  0127  200  FORMAT!///*  **  TFRMI NA  TED  WITH  ISTFP**,15,*  LASTEP=*,I5, 

1  *  TIME=*,1PF11.3,»  TMLAST**,E11.3,*  IFIN** ,13) 

SN  0)28  GO  TO  110 

SN  0129  320  STOP 

SN  0130  END 
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A.. 6.2  Subrouting  SfYVAP 

•  c  IM  IPl 

«C 


SN  0002 


i 

! 


SN  0003 


SN  000A 

SN  0005 

SN  0006 
SN  0007 


SN  0009 
SN  0010 
SN  0012 
SN  0013 

SN  0016 


SUBROUTINE  SR  TV  AH  I  OUT  I 

C***»*M  F.TAMANINl,  FACTORY  MUTUAL  RF SEARCH  CORE.  >  MARCH  1976  ••*•••*> 
C* 

C*  THIS  SUBROUTINE  COMPUTES  TEMPERATURE  AND  DENS  ITT  PROFILES  FOR 
C*  ONEDIMENSIONAL  UNSTEADY  HEAT  CONDUCTION  IN  A  SOLID  SLAB  OF 
C*  FINITE  THICKNESS  UNDERGOING  PYROLYSIS. 

C* 

C*  TWO  TYPES  OF  SURFACE  BOUNDARY  CONDITIONS  CAN  BE  HANDLED: 

C*  11  PRESCRIBED  TEMPERATURE  AT  IOP  NF AR 1  THE  SURFACE 

C*  1)  PRESCRIBED  SUREACE  HEAT  FLUX  INOT  INCLUDING  CONVECTIVE 

C*  HEAT  TRANSFER  OR  SURFACE  R ERADI AT IONI 

C* 

C*  THE  R ATE  OF  PYROLYSIS  IS  GIVEN  BY  A  FIRST  ORDER  ARRHENIUS  REACTIO 

C* 

C*  PYROLYSIS  GASES  DIFFUSING  THPOUC.H  THF  SOLID  ARE  ASSUMED  TO  8F 
C*  IN  PERFECT  THERMAL  CONTACT  WITH  THE  SOLID  AND  TO  BE  MOVING  IN 
C*  THE  0 1°  FC  T I ON  OF  OF CR FA  St  NG  SOLID  DENSITIES  OR  TO  MIGRATE 

C*  TOWARD  BOTH  FREE  SURFACES  IN  SUCH  A  MAY  THAT  THE  NET  PRFSSJRF 

C*  DROP  ACROSS  THE  SLAB  IS  7E«P. 

C* 

C*  SPECIFIC  HEATS  (ACTIVE  SPUD.  CHAR  AND  PYROLYSIS  GASES1  AND 
C*  THERMAL  CONDUCTIVITIES  (ACTIVE  SOHO  AND  CHARI  ARF  TREATED  AS 
C*  LINEAR  FUNCTIONS  OF  LOCAL  TEMPERATURE. 

C* 

C*  S.l.  UNITS  APE  USED  THROUGHOUT  (KG.M.SECI 

C* 

c*** ************************************************************ *•***« 

COMM0N/SPYCOMABC(2I  ,CPA,C«»C.CPG,CPAS,CPCS,CPGS,DARCY,DT  I  MAT, 

1  OTIME  t0T( M!N»0TSTF°»DX»EPS1 , EPS? , ERRMX ,HF LUX  1 , HFLUX?, 

?  HI ,H ?, IBC( 21 , IF  IN, I STEP,ITER,1TERID,ITERMX,KBC( 21, 

3  LA^TEP,MG(46),MG1 ,MG?,N,NP1 ,PF ,QPO ,01DT , Q20T , REL AX, 

4  RHQ(  44), RHOA ( 44 ) ,RHOF , RH OM , SU MI T ,T (44 1 , T B AR ( 44 1 , 

5  TBC( 2,1*0) ,TCA,TCC,TCAS,TCCS,TIME,TIMEAV,T IMEBC(15D|, 

6  TINF1,TINF2,TMLAST,TRNEG,T0K,XBC(2) 

DIMENSION  A( 441,6(44)  ,BCTR1 (21 , BCTR2 ( 2 ) ,BT ( 44| ,C ( 44  1 , 0(44), 

1  OBHO(  441 ,F ( 441 ,1  SURE (21 ,M(441 .SLOPE ( 2  I , T P( 44 1 , TSL CR  (  2 

2  TSTRK21  ,TSTR2(2l  ,U(44|  ,XOF  PTH  ( 2 1  ,Y0(  2 1 
REAL  M,hg,MG1,MS2 

C 

CHAPTER  000000  -  PRELIMINARIES  -  00003300000000 

C 

1  STEPR=  I  STEPMOUT 
IF(  ISTFP.GT.OI  GO  TO  310 
C 

CHAPTER  1  1  l - GEOMETRY,  CONTROL  INDEXES  AND  THERMODYNAMIC  VARIABLE 

C 

K«1 

IFtN.LT.44)  GO  TO  125 
MR  1 TE ( 6, 1  301  N 

130  EOPMATt  ^DIMENSION  OF  ARRAYS  IS  INSUFFICIENT  TO  HANOLE  •, 

1  13,'  GRID  POINTS*) 

I F IN* l 
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-l 

2.1  t 

,  J  N 

76 

)  ■aciotv  mutual  muo*  oowoiawn 

8*11.7 

SN 

0015 

•CTUiN 

SN 

0016 

125 

CONTINUE 

SN 

0017 

NP1*N«1 

SN 

0018 

ROEN-ALOG! l./ERRMXI 

SN 

0019 

0TIME-DT1M1N 

SN 

0020 

T  IME*TIME8C 111 

SN 

0021 

OX-OX/FLOATfN) 

c- 

- - 

-  INDEXES  FOR  CONTROL  OF  SURFACE  B.C.  <J*1,FR0NT;  x2,BACKI  — 

c 

KBC(JI«1,  TEMPERATURE  B.C.  ;  »2,  HEAT  FLUX  B.C. 

c 

IBCIJI:  FIRST  GRID  POINT  LEFT  OF  WHERE  TEMP. B.C.  IS  IMPOSED 

c 

IBCtJI~0  FOR  SURFACE  TEMP. B.C. 

SN 

0022 

00  110  J* 1 . 2 

SN 

0023 

X0EPTH1  J»«X8CU» 

SN 

0024 

IF(KBC< JI.F0.2)  GO  TO  105 

SN 

0O?6 

IFIXBC(  Jl.EO.O.)  GO  TO  106 

SN 

0028 

XBC  <  J )* XBCI  J ) /D X 

SN 

O02<» 

TBC<J1*1MNT(XBC(JI1 

SN 

00  30 

XBC(JlxXBCIJI-AINT<  XBC  I4»> 

SN 

0031 

lFIJ.EO.il  GO  TO  110 

SN 

0033 

IBC ( J 1XNP l-I9C< Jl 

SN 

0034 

xbc(  j i* i.-xe:< ji 

SN 

0035 

GO  T(1  110 

SN 

0036 

105 

xbci  Ji-o. 

SN 

0037 

106 

IBC( Jl-0 

SN 

0038 

110 

CONTINUE 

L* 

h  in  1  nun  r  T  imjl  T  ii  5  i  mrcnM  i  un  c  —  —  —  —  — 

SN 

0039 

CPW=CPA 

SN 

0040 

TCWXTCA 

SN 

O041 

RHOFM1=1.-OHOF 

SN 

0"42 

TL*  l.E  30 

SN 

0063 

1FIPF.E0.0.I  GO  TO  115 

S*' 

0065 

TRNFG«-TRNEG/8314. 

S*I 

0066 

P  F=PF /RHOFM 1 

SN 

0067 

Tl*-TRNEG/l ALOGI 100.6PF*CPW*RHOW/TCW1»2. *ALOG(DX*Fl OAT  INI  1 1 

SN 

0068 

MR  I TEI  6»  1 20 1  TL 

SN 

0069 

120 

FORMAT! /•  PYROLVSI  S  CALCULATION  IS  NOT  PERFORMEO  FOR  TEMPf 

1TURES  LESS  THAN' «F  B. 2v  *  OEGK ' /75 X , •*** •// 1 

SN 

0050 

115 

C0NT1NUF 

V 

jr  tL  X  i  XL  nr  A  1  j  A<iU  intKWAL  L  U'NUUt  I  XVII  1  to  —  —  —  —  —  —  —  —  —  —  —  — —  — 

SN 

0051 

THD  I FF*  TC M/RHOM/C  PM 

SN 

0052 

RCDX=RHOW*CPW*DX 

SN 

0053 

CPA=CPA/CPM 

SN 

00  56 

CPAS=CPAS/CPM 

SN 

0055 

CPC=CPC/CPM 

SN 

0056 

CPC S*CPCS/CPW 

SN 

00  57 

TCMOHX* . 5*TC  M/D  X 

SN 

0058 

TCA*TCA  /TCW 

SN 

0059 

TCAS*TCAS/TCM 

SN 

0060 

TCC»TCC/TCM 

SN 

0061 

TCC  S*TCCS/TCM 

SN 

0062 

CPG04«.25*CPG 

SN 

0063 

IF!CPG$*CPG.NE.O.i  CPGS*.5*CPGS/CPG 

SN 

0065 

EPS1*FPS1*5. 669E-8 

SN 

0066 

EPS2«EPS2*5.669E-6 

94 


SN  0067 
SN  0068 
SN  0069 
SN  0"70 
SN  0"71 

SN  OP 73 
SN  0075 
SN  0076 
SN  0P78 
SN  0079 
SN  0080 
SN  0081 
SN  ''PS? 
SN  0083 
SN  0"84 
SN  0085 
SN  0086 
SN  0087 
SN  0088 
SN  0089 
SN  0090 
SN  0091 
SN  0092 
SN  0093 
SN  0094 
SN  0095 
SN  0"96 
SN  0097 
SN  0098 
SN  0099 


SN  0100 

SN  0101 
SN  01"2 
SN  0103 
SN  0105 
SN  0106 
SN  0107 
SN  0108 
SN  0109 
SN  0110 
SN  0112 
SN  0113 
SN  0114 
SN  0115 
SN  0116 
SN  0117 
•SN  OH® 


SN  0119 

i 


FACTOAY  MUTUAL  KStAICM  COVOtATIOM 

X 1011.7 

C 

CHAFffK  *22—  INITIALIZATION  OF  AKA  AYS  AM)  OTlCK  VAN  I  ABIES  ?  ? 

IF IN-0 

00  210  J=l,2 
TSLCR!  Jl-l. 

SLOPE! J 1= .01 

IF! TBCI J, 1).NE.T8C!J,211  SLOPE! Jl-ITBC! J,21-TBCI J,1 ) 1/ 

1  ITIMEBCm-TINEBCmi 

IF! XBC! Jl .FO.O. I  GO  70  210 
X0=  .5*XDE  PTH!  Jl  /SQP T! THOI F F *DT I  MAX  1 
IF!  XO.GT.  10.  1  XD  =  10. 

SLOPE!  J  I* SLOPE!  Jl  /!  1 . ♦  2.  * XD*XD1  /EBFC I XD 1 
210  CONTINUE 
I  SUP  F I  1 1  =  1 
1  SUP  F I 21*NP 1 
0  1DT=0. 

0  2DT=0  • 

00 1DT=0 . 

002DT*0. 

MG1=0. 

MG2=0. 

C  1=0. 

C  2=0  . 

00  220  1=1, NP1 

T|  I  1=T0K 
TP!  I  l  =  T"K 
8  HO  (  I  1*1. 

RHOA!  I  1*1. 

OP HO! I  1  =  0. 

MG!  I  1  =  0. 

C(  I  1  =  0. 

01  I  1  =  0. 

22"  CONTINUE 

CHAPTER  333333  - EVALUATION  OF  BOUNDARY  VALUES -  3333333 

C - LINEAR  INTERPOLATION  OF  BOUNDARY  CONDITIONS  - 

31"  CONTINUE 

T!"EAV=TIME*. 5*0TI*E 
T|MF0T=TIME*0TINE 

IF!  TIMEAV.LE.TIME8CIKMI  1  GO  TO  311 
K=K  ♦! 

GO  TO  31" 

311  CONTINUE 

00  315  J=l,2 
BC! JI*T1ME0T 

1FIKBC! JI.E0.2I  BCIJ1*T!NFAV 

BC! Jl*!  BC!  JI-TINEBC! K1 l/t Tl MFBC I K*1 1 -T I MEBC! XI 1 
315  BC! J)*T8C!JtKI*8C! Jl*! TBCI J,K»1 l-TBCf J,KII 

IF  (KBC(l).EQ.l)  T?(l)  =  BC ( 1) 

IF  (K3C(2).EQ.l)  TP(NPl)  =  BC(2> 

HFLuxi-ncm 
HFLUX2«BC!21 
I TER*0 


C 

C - CALCULATES  SURFACE  TEMPERATURES  I  FRONT ,  J*1 ;  BACK,  J*2I 

C  WHEN  TEMPERATURE  NEAR  IAN0  NOT  ATI  THE  SURFACE  IS  GIVEN 

C  AS  BOUNDARY  CONDITION 


325  CONTINUE 
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SN 

01 

20 

SN 

0122 

SN 

01 

23 

SN 

01 

25 

SN 

01 

27 

SN 

01 

28 

SN 

01 

29 

SN 

01 

30 

SN 

0132 

SN 

0133 

SN 

01*4 

SN 

"136 

SN 

0138 

SN 

0140 

SN 

0141 

SN 

0142 

SN 

0144 

SN 

0146 

SN 

0148 

SN 

0150 

SN 

0152 

SN 

0153 

SN 

0155 

SN 

0156 

sv 

0158 

SN 

0159 

SN 

0160 

SN 

0161 

SN 

0162 

SN 

0164 

SN 

0166 

SN 

0167 

SN 

0168 

SN 

0169 

SN 

0170 

SN 

0171 

SN 

0173 

SN 

0174 

SN 

0176 

SN 

0177 

SN 

0179 

SN 

0180 

SN 

0181 

SN 

0182 

SN 

0184 

SN 

0186 

SN 

0187 

SN 

0189 

SN 

0191 

IFflftC<tl«f4.0.AM0.I*Cm.E0.0l  60  TO  220 
00  350  J-1.2 

IFI ISC! Jt.EO.O)  GO  TO  350 
I FI  ITER .GT. 01  GO  TO  33  0 
TSTR1IJI-TI  I  SURF  I  Jl I ♦SLOPE!  JI*DTIME 
TPI ISURFI Jl J  =  TSTRl I J) 

GO  TO  350 

330  IFI  ITER.GT.il  GO  TO  340 
C 

C - FIRST  ITERATION 

BCTP1I  J)“TPI  IBCI  JIIaXBCI  Jl*|Tt>|  IBCIJlAll-TPIIBCI  Jll  1 
XO*XOEPTHI Jl*. 5/SQPTI THDI FF *OTI ME  I 
IFI  XO.GT.  10.1  X0-10. 

IFI XO.LE. 1. 51  FACT“1.-.6402*XD 
IFI XO.GT  . 1.51  FACT=. 5642*F XPI -XD*XD> /X 0 
FACT*|  l.*2.*XD*XDl  *FACT 
YOI J)*FACT 

IFIFACT.IT..01I  F  AC  T  =  .  01 
IFI ISTEP.EQ.OI  BCTR1I J1“BCI Jl 

IFIBCTRUJI.NE.BCUn  TSLCR  I  J)  =  SLOPF I  J)  ♦  I  BCI  J I-BCT  PI  IJ 
1  /OTIMC/FACl 

IFIBCTR1I Jl.tQ.BCI Jll  TSLCRIJ1“1.1*SL0PE!JI 
IFI SL  IMI T.GT. 0. I  GO  TO  331 

YDI J1*-SL IMI T*0TSTEP*DX*FL0AT|N1*. 5/SQPTI THDI FF*OT I  ME  1 
IFI TSLCR! Jl /SLOPE  I JI.LT. I.)  YDIJI=10.*YDIJ1 
FACT*ABSI TSLCR I Jl- SLOPE  I Jll *DTIME/VDI Jl 

I  F|  FACT .GT .  1.  I  TSLCR!  Jl  “SLOPE  I  JIM  TSLCR  IJ1 -SLOPE!  Jl  1/TACT 
GO  TO  332 

331  CONTINUE 

YOI JI*SL!MIT/ll.-.1*AL0GIVDIJl  II 
FACT=TSLCRI Jl/SLOPEI Jl 

IFI FACT .L  T. 0. . ANO.-F  AC  T.GT. YO I J 1 1  TSl CR  l  Jl  =-Y 01  J) *S  LOPE  I d  l 
IFI FACT.GE .0. . AND. AB SI  FAC T-I .  I  .  GT. VOC Jl  I 
1  TSLCR1JI* SLOPE  I J1*I1.*SIGNIYDI Jl ,FACT-i.l! 

332  CONTINUE 

T  STB  21  J  1=  Tl  ISUPFI  Jl  IMSlCRI  Jl  *DTI  ME 
TPI  ISURFI  JII-TSTR2IJI 
GO  TO  350 
C 

C -  S ECONO  ITERATION 

340  BCTR2IJIxTP(IBCIJll*XBClJl*ITPIlBCUm  I  -T  P 1 1  BCI  J 1 1  I 
IFIBCTR2I  JI.NE.BCTR1IJ11  TPI  I  SURFIJl  1*TSTR2I  JIMTSTR1IJ  I-TSTR2I. 

1  *IBCI JI-BCTR2IJ1 l/IBCTRl «  Jl -BCTR2 I Jll 

TSLCR I J 1*1 TPI I SURFIJl 1 -TIT  SURF! Jl  1 1/OTI  ME 
IFI SLIMIT.GT.O. I  GO  TO  341 
FAC  T-ABSI  TSLCR I Jl -SLOPE  I Jl  1  *OTI  ME /YOI J» 

I  F(  FAC  T.GT.  1.  I  TSLCR  I  Jl*  SLOPE  IJIMTSI  CPI  Jl-SLOPEIJ  1 1/FACT 
GO  TO  342 

341  CONTINUE 
FACT-TSLCRI Jl /SLOPE! Jl 

I  FI  FACT. LT. 0.. ANO.-F AC T.GT. YOI Jl I  TSLCR I Jl *-YDI Jl *S LOPE l J 1 
IFIFACT.GE.C. .ANO.ABSIFACT-l.l.GT.VOI J) 1 
1  TSLCR! J) “SLOPE 1J)*II.*SIGN|Y0IJ1 , FACT- 1.1) 

342  CONTINUE 

IFI 1STEP.GT .0)  SLOPE  I J) “TSLCR I Jl 
IFI  SLOPE!  Jt.EO.O.  I  SL0PEIJI*.91 

IFIABSI SLOPE  I J 1 1.LT..01I  SLOPE! Jl “SIGN!. 01 .SLOPE  I Jl  1 


96 


SN 

0193 

SN 

0194 

SN 

0195 

SN 

0196 

SN 

0197 

SN 

0198 

SN 

0200 

SN 

0  701 

SN 

0202 

SN 

0?o<, 

SN 

0  ?0  6 

SN 

O?07 

SN 

0209 

SN 

0?1O 
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f»l  ISUIFUIt-T(ISUM(JI1»k»E(J>tOTIN( 

350  COMYnWE 
320  CONTINUE 

CHAPTER  44444444 - BEGINNING  OF  LOOP - 

C -  COMPUTES  DENSITY  INCREMENTS  - 

400  CONTINUE 

OP  410  1=1, NP1 

I  FI  ITER  .E0.0)  RH0A(  n=RH0II  )».S*DRH0m 
T  BAR III=.5*ITP(I|*T(I|| 

PR HOI  I 1  =  0. 

IFITBARI I I.LE.TL)  GO  TO  405 
1FIPH0I  I  I.EO.RHOF)  GO  TO  405 

DP  HOI  I l=DTIME*PF*I RHOF-RHOAII » I *E X P (TRNEG/T B AR (II ) 
IFIRHOI  ll«-OPHOIII.LT.RHOFI  ORHO  (  I  I  =PHOF-RHOm 
405  RH5AI  I»  =  RHO(  II*.5*0RH0m 
410  CONTINUF 


4444444444 


-COMPUTES  NEW  DISTRIBUTION  OF  GASEOUS  FLUX,  ASSUMING  That 
FLUX  IS  IN  THE  DIRECTION  OF  DECREASING  DF  NS  IT  V  ( DARCY .N e. 1 .  1 
OR  THAT  IT  FOLLOWS  DARCY'S  LAW  IDARCY. EQ. 1. » 


SN 

021  1 

IFIRHOAI 1I.FQ.1..AND . RHOA (NP1I.E0.1.)  GO 

TO  435 

SN 

0213 

J*NP  1 

SN 

0214 

MG(NPH  =  0. 

SN 

0215 

N  G2=Q . 

SN 

0216 

DO  4  70  I  =  2, NP 1 

SN 

0217 

ICPL=NP  !♦  1-  I 

SN 

0218 

MGI  ICPL)  =  MG<  ICPLMl-DRHO<!CPL*n*DX*RHOW/DTIMF 

SN 

o  ?  1  9 

I  F(  ICPL.EO.NI  MGI ICPL1=.S*MG( ICPLl 

SN 

"221 

IFIRHOI  ICPL  I.GT.RHOAIJII  J=!CPL 

SN 

0223 

470 

CONTINUF 

SN 

"274 

M  G1  =  MGI  1 1-DRHOI 1 ) * . 5*0  X*°H0w/DT I  ME 

SN 

0775 

IFIOARCY.NE.l.)  GO  TO  415 

SN 

0227 

MG2= .25*1 MGI- MGI 11 -MGI NI *MG I NP1 ) I 

SN 

0778 

DO  416  1=1, N 

SN 

0779 

416 

MG7=MG2*MGI  I  1 

SN 

0730 

MG2=MG2/FLOAT(N) 

SN 

0231 

GO  TO  417 

SN 

0732 

415 

CONTINUF 

SN 

0233 

IFIJ.EO.NPII  GO  TO  435 

SN 

0735 

IF! J.NE.l J  MG 2=. 5*1 MGI JI*NG(J-l 1 I 

SN 

0737 

IFIJ.EO.il  MG 2* MG  1 

SN 

0239 

417 

CONTINUE 

SN 

0740 

00  430  1=1, NP1 

SN 

0741 

430 

MGI  I  )  =  MGI  I  1  -MG 2 

SN 

0742 

MG1  =  MG 1-MG  2 

SN 

0243 

MG?=-MG(  NP1 1 

SN 

0244 

435 

CONTINUE 

L 

C- 

..... 

-  COMPUTFS  AVERAGE  VALUES  OF  SPECIFIC 

HEAT  - 

SN 

0745 

00  440  1=1, NP1 

SN 

0246 

A 1  I  1*PC  DX/DTT  ME 

SN 

0747 

IFIRHOAI! I. EQ.l.. AND. CPA S. E 0.  0.)  GO  TO  440 

SN 

0249 

CORR=CPA*CPAS*I  TBARII!  -T0K1 

SN 

0250 

IFIRHOAI 1 I.EO.l.)  GO  TO  441 

SN 

0252 

CORR*I I RHOA 1 I l-RHOFI *CORR*RHOF*Il .-RHOAII) l*ICPC*CPCS* 

1  1 TBARII l-TOKII | /RHOFMl 

SN 

0253 

441 

All  l-CORR*A( I | 
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SN 

0254 

440 

CONTINUE 

SN 

0255 

MlKHIIll 

SN 

0256 

A(NPI)*.«.*A(NP1 1 

L 

c- 

- COMPUTES  AVERAGE  VALUES  OF  THERMAL  CONDUCTIVITY - 

SN 

0257 

DO  445  I»1,NP1 

SN 

0258 

81  I  »*TCMDHX 

SN 

0250 

IFIRHOAU  I.EQ.1..AN0.TCAS.E0.0.I  GO  TO  445 

SN 

0261 

CORR  =  TC  A*TC  AS*(  TBARI I l-TOK) 

SN 

0262 

IFIRHOAU). FQ.l.)  GO  TO  446 

SN 

0264 

CORR  *(  ( RHOA  (  I  )-RHOF)*CORRMl.-RHOACm*tTCC*TCCS*|TBARm- 

1  T0KDI/RH0FM1 

SN 

0265 

446 

8(1 )=CPRR*8(l  I 

SN 

0266 

445 

CONTINUE 

SN 

0267 

00  4  50  1*1, N 

SN 

0260 

8(1)*. 5*181 11*81141(1 

SN 

''260 

450 

CONTINUE 

l 

SN 

0270 

IFIMG1.E0.0..AN0.MG2.EQ.0.I  GO  TO  463 

SN 

0272 

00  455  1*1, N 

SN 

0273 

Cl  I  >  =  Hr.(  1  )*CPG04 

SN 

0274 

IFICPGS.EO.O.)  GO  TO  455 

SN 

0276 

C  I  I  l  =  C( I  )*( 1 .  ♦( 1 TBARI 1 )♦ TBARI !♦ 1 ) >  *. 5-T 3K) *CPGS 1 

SN 

0277 

455 

CONTINUE 

SN 

0278 

C 1«MG1*CPG0  4*( 1 .♦( TRAP  1 1 l-TOK) *CPGS> 

SN 

0279 

C2*-MG2*CPGD4*( 1 .♦( TBARI NP1) -TOK) *CPGS) 

SN 

0280 

460 

CONTINUE 

L 

c- 

-  COMPUTES  ENERGY  SOURCE  DUE  TO  PYROLYSIS  IN  THE  SOLID  - 

SN 

0281 

IFIPF.EO.O.)  GO  TO  470 

SN 

02«3 

00  465  1=  1 , NP 1 

SN 

0284 

01  I  (  =  0. 

SN 

0285 

IFIORHOI  D.EO.O.)  GO  TO  465 

SN 

0207 

Qp=TRAR I T l-TOK 

SN 

0288 

0P  =  QP*(CPA*.5*C  PA  S*QP-RHOE*(  CPC*-.  E  *C  PCS*QP) ) 

SN 

02P9 

OP  =  QPO*CPW*QP/R  HOF  Ml 

SN 

0290 

01  I  )  =  -QP*DRHO( I |*D X*RHOW/OTIME 

SN 

0291 

465 

CONTINUE 

SN 

0292 

D(  1  )=. 5*0(1) 

SN 

0293 

D I  NP  1  )* ,5*DI NP1 ) 

SN 

0294 

470 

CONTINUE 

CHAPTER  555555  - TRIDI AGONAL  MATRIX  OPERATIONS -  5555555 

SN 

0295 

00  510  1*1, N 

SN 

0296 

F(  I )*-B 1 1 l-CI I) 

SN 

0297 

510 

CONTINUE 

SN 

0298 

U 1 1  »  =  B I  1)-C(1)*2.*C1*.5*(H1*EPS1*TBAR(1  )**3) 

SN 

0299 

DO  515  1*2, N 

SN 

0  300 

Ul  I  1*81  I  1*81 I-1I-C(I)*C(I-1) 

SN 

0301 

515 

CONTINUE 

SN 

0302 

UINP  1)*B(N)*:(N)-2.*C2*.5*(H?*EPS2*T8AR(NP1  )**3) 

SN 

0303 

00  520  I-2.NP1 

SN 

0304 

Ml  I  )-C( I-l)-B(I-I) 

SN 

0305 

520 

CONTINUE 

SN 

0306 

BTIll-tAI  l)-U( l))*T(l)-Fll)*T(2)*4.*T3K*(Cl-C(l) )*0(1)*H1*TINF1 
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SN 

0307 

00  f»  1*2*11 

SN 

0308 

1  4.*T0K*ICII- , -Clll  )»D(I) 

SN 

0300 

525 

CONTINUE 

SN 

0310 

BTCNPll  —  MCNP11  A<  NPl  1  -U(  NP1 1  1  *T  <  NP1  »*4.*T0K*  1 C  1 N  >-C? )  ♦ 

1  DCNP1»4H2*TIHF2 

SN 

0311 

00  530  1*1, NPl 

SN 

0312 

UC  1  l*UC  !  * ♦*  1 1  » 

SN 

0313 

530 

CONTINUE 

L 

C- 

... 

-  MAKES  APPROPRIATE  CHANGES  !N  COEFFICIENTS  WHEN 

A  SURFACE 

c 

TEMPERATURE  ROUNOARY  CONDITION  IS  GIVEN 

c- 

... 

FRONT  SURFACE 

SN 

0314 

!E<  KBC(  ll.EO.il  GO  TO  540 

SN 

0316 

8  T(  1 1*8  T(  1 1 ♦BC  C  11 

SN 

031  7 

GO  TO  545 

SN 

0318 

540 

CONTINUE 

SN 

0310 

BT(  1  1=8 T(  n-Ulll*TPUI 

SN 

0320 

BTC  2»*BTC  71-M(  21*TPm4RTU> 

SN 

0371 

U1  2 1  =  UI  2)*F  C  1  1 

SN 

0377 

U(  11—1. 

SN 

0373 

MI71—I. 

SN 

0374 

545 

CONTINUE 

c- 

— 

BACK  SURF  AC  F 

SN 

0325 

IFCKBCC 2l.E0.il  GO  TO  550 

SN 

0327 

BT(NPll*BT|NPl)4BCm 

SN 

0378 

GO  TO  555 

SN 

03  ?o 

550 

CONTINUE 

SN 

0330 

ft  T ( NPl 1=8 TI NPl 1-U1NP1 1 *TP1NP1  1 

sv 

0331 

6  T(N )*PT1N1-F(N»*TP( NP11»9TINP1  > 

SN 

0332 

U<NI*U(N|*PMNP1I 

SN 

0333 

U(NP 1 1*-  1  • 

SN 

0^4 

FCN1*-1. 

SN 

0335 

555 

CONTINUF 

L 

c- 

— 

- SOLVFS  THE  TR 1 OT  AGONAl  MATRIX  AND  PUTS  NEW  T F MPER  AT ORES 

c 

IN  ARRAY  ftTU  » 

SN 

0336 

00  560  1*2, NPl 

SN 

0337 

M(  I 1*M(  I  »/U( I-l  ) 

SN 

0338 

560 

U1 I 1  =  UC  I 1-N1 I l*F( 1-11 

SN 

0*39 

DO  565  1*2, NPl 

SN 

0340 

565 

BT(  I1*BT(  Il-MU  1  *BT<  I  -1 1 

SN 

0341 

BT(NP1  i-btcnpii  /UC  npu 

SN 

034? 

00  570  1*1, N 

SN 

0343 

J*NP 1- I 

SN 

0344 

570 

BTC  JI*CBTIJ1-FC  J)*BTC  JM11/UCJ) 

c- 

— 

-  MAKFS  APPROPRIATE  CHANGES  IN  SURFACE  VARIABLES 

WHFN 

c 

SURFACE  BOUNDARY  CONDITION  IS  ON  TEMPERATURE 

SN 

0345 

IFCKBCC  11.E0.2I  GO  TO  585 

i 

SN 

0347 

HFLUX1«BTUJ 

SN 

0348 

BTC  II-TPC  11 

! 

SN 

0349 

585 

CONTINUE 

t 

SN 

0  350 

IFCKBCC 21. EO. 21  GO  TO  590 

1 

SN 

0352 

HFLUX2*BT(  NPl  1 

1 
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$N  0353 
$N  0354 


SN  0355 
$N  035e> 
SN  0357 
SN  0358 
SN  0360 
SN  0362 
SN  0063 
SN  O06<, 
SN  0365 
SN  0066 
SN  0367 
SN  0068 
SN  O06<» 
SN  0370 


^ N  0371 
SN  0073 
SN  0375 


SN  0376 
SN  0378 

SN  0080 
SN  03fl? 


SN  0384 
SN  0386 
SN  0387 
SN  0380 
SN  0390 
SN  0391 
SN  "392 
SN  0393 
SN  0394 
SN  0395 


SN  0396 
SN  0397 
SN  0398 
SN  0399 

SN  0400 
SN  0401 
SN  0403 
SN  0404 


8T|*Pil>T»tNMt 

590  CONTINUE 
C 

CHAPTER  666  - CONTROL  Of  NUNBER  OF  ITERATIONS  ANO  STEP  SI7F -  6  6  < 


C 

C -  CHECKS  FOR  MAXIMUM  ERROR  BETWEEN  SUCCESSIVE  ITERATIONS 

C  ANO  UPDATES  TPIII  ARRAY 


ERP«0. 

CHG=0. 

00  610  I-1.NP1 

IF!  I.E0.1.0P.I.E0.NP1I  GO  TO  607 
IFIABSIBTI Il-TPI Ill.LT.ERR)  GO  TO  607 
ERR*ABSI BTI ll-TPl I  I  I 
IER=  I 

607  CONTINUE 

BETAN1  =  TP( I  I 
TP(  l  >  =  BT<  I  » 

BTI  I  )=RETAM1 

C  HG=AM A  XlICHGtABSI  TP(I)-T|I  1)1 
610  CONTINUE 

I  TER* ITER  ♦  ! 

C 

C - 00 ES  3  PRELININARY  ITERATIONS  WHEN  B.C.  IS  NOT  AT  THE  SURFA 

IF  ( IBCI 1 1.FQ.O.ANO.IBC12I.EO.OI  GO  TO  605 
IF( ITER. IT. 31  GO  TO  325 
605  CONTINUE 
C 

C - CHECKS  FOR  REQUIREO  ACCURACY  ( ERR NX  I  - 

!F(  ISTEP.EO.OI  GO  TO  615 
IF  (PF  +  CPAS  +  TCAS  +  EPS1  +  EPS2  .EQ.  0. 

1  .AND.  IBC(l)  +  IBC(2).EQ.  0)  GO  TO  615 

IFIERR.LT.FRRNX.0R.ERR.LT.ERRMX4CHG)  GO  TO  615 
IF!  ITER .EQ.ITERNXI  GO  TO  625 
C 

C -  RELAXATION  OF  ESTIMATED  TEMPERATURE  PROFILE  - 

IFIRFLAX.EO.O.I  GO  TO  400 

BFTA=ERR/FBRMX 

IFICHG.GT.1.1  BETA-BETA/CHG 

BFTA  =  AM  INK  1 .  ,  A  LOG  I  BE  T  A)  /R0EN1 

BETA=RELAX*BETA 

BFTA*  1* 1,-BET  A 

00  640  !*1'NP1 

TPI  n=BETA*BTII  )*BFTAMl*TP(n 
640  CONTINUF 
GO  TO  400 


C 

C -  TIME  STEP  IS  HALVED  WHEN  TOO  MANY  IITERMX1  ITERATIONS  ARE 

C  REQUIREO.  EXECUTION  STOPS  WHEN  OT I  ME  IS  LESS  THAN  OTIMIN 


625  CONTINUE 

ERR  *ERR /CHG 

WRITE (6, 6111  ITER, ISTEP,IER, ERR 

611  FORMAT!  •  **•  MORE  THAN*, 13,  •  ITERATIONS  REQUIRED  AT  ISTEP*',I4, 
1  •  /  TEMP.  RELATIVE  ERROR  AT  I-', 13,'  WAS  •.IPE11.41 
0TImf«.5*DTIME 

IF! OTIME.GT. OTIMIN)  GO  TO  310 

IFIN*2 

GO  TO  630 
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SN  0406 
SN  0407 
SN  0408 
SN  0409 
SN  0410 

SN  0411 

SN  0412 
SN  *413 
SN  0414 
SN  0415 


SN  0416 
SN  0417 
SN  f>410 
SN  0419 
$N  0420 
SN  0421 
SN  r>4?2 
SN  0423 
SN  0424 
SN  0426 
SN  0427 
SN  0429 
SN  0430 
SN  0432 
SN  0433 
SN  0434 
SN  0436 
SN  0437 
SN  0438 
SN  0440 
SN  0441 
SN  0442 
SN  0443 
SN  0444 
SN  0445 
SN  0446 
SN  0447 
SN  0448 
SN  0450 


SN  0451 
SN  0452 
SN  0454 
SN  0456 
SN  0458 
SN  0459 


619  CONTINUE 


C - COMPUTES  TOTU  ENERGY  WHICH  ENTERED  THE  SLAB  THROUGH 

C  FRONT  (01DT1  ANO  BACK  (Q2DT)  SURFACE 


Q1DT*Q1DT*0Q1DT 

Q20T«02DT*DQ2DT 

TBAR ( 1)*.5*(TP( 1 !♦ T€ 1 » 1 

T8AR ( NR  1 I*. 5*1 TPI NPl  1  ♦ Tl  N*1 1 1 

DQ1DT*HFLUX !♦«!*( TINE l -TBAR (l I » -EPSl *TBAR(1» **4 
1  -4.*C1*ITBAR|1»-T0KI 

DQ2DT*HFLUX 2*H2*I TINF2-TBAR(NP1>I-EPS2*TBAR(NP1I**4 
1  ♦ 4. *C2*( TBAR( NPl I -TOKI 

DQ10T*.5*0Tt  HE *001 OT 
D02DT* . 5*0T 1  HE  *0Q2DT 
0 1DT*Q10T*DQ10T 
020T*020T40020T 


C 

C - UPDATES  Tl  HE  AND  DENSITY.  CALCULATES  NEW  STEP  SUE  AND 

C  MAGNITUDE  OF  TPdl  AND  DRHOU)  AT  NEXT  STEP 


T1ME*TIME0T 

C  Kj=0  . 

DP  645  1*2, N 
FRP=A8S( TP( 1 1-TII I  1 
CHG=CHG^ERP 
645  CONTINUE 

C  HG=CHG/F  LOAT  (  N-l  1 
EPR=DTSTEP/CHG 
I F ( ERR .GT .1.1  ERR=1./ERR 
ERR  =  1  .-ERR 

TF(  ITER .GT.ITERID1  ERR *ERP*FL OAT (I  TER  I D I /F LOAT I ITER1 
OT 1*6  =  A  MI Nl( DTI  MAX, DTI  ME *( DTSTE P/CHG1  **ERR1 
I  Ft  ITFR.LE.ITER10I  GO  TP  635 
DTIHE*DTIME*SQPT(FLOAT(I TERID1 /FLOAT ( IT  ER1 ) 

635  CONTINUE 

TF< DTIHE .GT.DTIHINI  GO  TO  630 
IF IN=  3 

630  CONTINUE 

IF(  TIHE40TIHF .GT.THLAST1  DTI  HE* T*l AST-T I  HE 
CHG=DT!HE/2./(TIHE-TIHEAVl 

DP  620  I*1,N»]"*' - - - —IF  (ISTEP.EQ.O)  CHG  -  0. 

TBAR(  l  )=.5*(  TP(  1IHIII  I 
ERR  =T( I  1 
T (  I  l  =  TP( I  I 

TPI  I)*T(!  )*CHG*( Till -ERR I 
R HO ( I )*RHO(  1 1 *DRHO( 1 1 
DRHPI I 1*CHG*DRH0( I ) 

I F ( OP  HD ( I 14PH0I I I.IT.RHOF)  DRHPI T) *RHPF-RHO(I> 

620  CONTINUE 
C 

CHAPTER  777777  -  END  OF  LOPP  -  7777777777777 

C 

C -  CHECKS  WHETHER  TO  RETURN  TO  MAIN  FOR  OUTPUT  - 

ISTEP-ISTEPM 

I F( TIME.GE.THLAST.PR.I STEP.EO.lASTEPl  RETURN 

IFI ISTEP.F0.ISTEPR1  RETURN 

IF!  IF  IN .NE .01  RETURN 

GO  TO  310 

ENO 
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SN  0002  SUBROUTINE  OUTRUTIIMI  WT| 

SN  0003  COMNON/SPVCOM/BC!  2)  , CPA, CPC  • CPG,CPAS,CPCS,CPGS. DARCY,  DTIMAX, 

1  DTIME.DTIMIN,0TSTEP,DX,EPS1  ,EPS2,ERRKX,HFLUX1,HFlUX2, 

2  HI, H  2.  IBC  12) ,IFIN,ISTEP,ITER,ITER!D,ITERMX,KBC!2), 

3  LASTEP »MG! 44) ,MG1 ,MG2,N»NP1 , PF ,QP0 , 01 DT , 02 DT , PEL  AX, 

4  RH0(44),RH0A(44t ,RHOF  ,  RHOW.SLI  MI  T  ,T  144  )  ,  T  BAR  144 )  , 

3  T8CI  2,  150)  ,TCA  , TCC , TC A S, TCCS ,Tl ME ,T l ME AV ,T  IMEBC! 150), 

6  TINF 1,T1NF2,TMLAST,TRNEG«T0K,XBC!2) 


SN  0004  REAL  NG,MG1,MG2 

tSN  oo05  REAL'S  LAS( 3) /4HTEMP .7H0ENSI TV, 8HMASSFLUX/ 

C -  COMPUTES  FRACTION  OF  INITIAL  WEIGHT 

SN  O""#,  R  KJBAR* 1. 

SN  0007  I F( R HO A ( 1I.EQ.1..AND. RH0A1NPI I.EQ.l.)  GO  TO  20 

SN  0009  RH0BAR«.54PH0A!  1  ) 

SN  0010  DO  21  1-2, N 

SN  oo  1  1  R  HT1BAR*RH08ARaRH0A!  I  ) 

SN  0012  21  CONTINUE 

SN  0013  RHOBAR-RHOBAR*. 5*RH0A(NPt) 

SN  0014  R  HOBAR»RHOBAR /F  LOA  T( N) 

SN  0015  20  CONTINUE 

C -  COMPUTES  HEAT  FLUXES 

SN  0016  QC0NV1-H1*!  TINF1-T8ARU)) 

SN  0017  RERA01«EPS1*T8AR!  11**4 

SN  00 1 R  QC0NV2«H2*< TINF2-TBAR1NP1I) 

SN  0019  RERA02«EPS2*TBAR!NP1I**4 

C -  PRINTOUT  OF  STATION  VARIABLES 

SN  "020  WP I TE! 6,  102 )  ISTEP.KBCUl  ,KBC1?)  ,ITEP,DTTME 

SN  0021  102  FORMAT!//*  ISTFP=*,I5,*  KBC 1 1  )  *  *  1 12  ,  *  KBC  1 2 1  »•  ,  I  2.  *  ITER**, 12, 

1  •  0T|ME-*»1PE10.3I 

SN  0022  WPTTE16.103)  MG  1 , MG2 , RHPBAR 

SN  0023  103  FORMAT!*  MG1 * • , 1  PE  1 0. 3 , •  MG2* • , 51 0. 3 , *  RHOBAR  =  • , 0 PF8 .5  I 

SN  0024  WRITEI6.104I  TI ME A V, HF LUX1 , OCON VI ,RFR A01 ,01 OT , 

1  HFLUX2,0CPNV2,RERAD2,020T 

SN  "025  104  FORMAT! •  Tl ME  A V-* ,F 8. 2 , *  OR AOl * • , 1  PEI  0. 3 , •  QCPNV 1 *• , f 10 . 3, 

1  •  RE®  AO  1*  * ,E1 0. 3, •  010T«*,F10.3/16X,*  QRAD2** ,E10 .3, 

2  «  QC0NV2*  *,E10.3,*  PERAD2** ,E10.3 , •  Q2DT** ,F10.3) 

c -  PRINTOUT  OF  BOUNDARY  CONDITION  MATCH 

SN  0026  IF! XBC! 1) .EO.O.I  GO  TO  23 

SN  0028  TSCC»T1  IBC!  11 )♦ XBC ! 1 1 *! T ! I BC! 1 ) *1 ) -Tl IBCC 1 1 > I 

SN  0"29  WRITE! 6, 105)  BC11I, TBCC 

5N  003O  105  FORMAT! •  FRONT  B.C . 5 • ,26X,FB. 3,7X,F8.3) 

SN  0031  WP I TE1 6, 106) 

SN  0032  106  FORMAT!  •♦*, 15X, ‘TEMPERATURE  SHOULD  BE : •  ,1 2X , • IS : •  ) 

SN  0033  23  CONTINUE 

SN  0034  IF! XBC! 2I.E0.0. )  GO  TO  24 

SN  0036  TBCC«T! IBC! 2) )♦ XBC !2J*!T!IBC!2)*1)-T!IBC!2))I 

ISN  0037  WR1TF16.107)  602), TBCC 

’  SN  0038  107  FORMAT!*  REAR  B.C . : • , 2 7X,FB, 3 , 7X ,F8 . 3 ) 

■’SN  0039  IF1X8CI  D.EQ.O.)  WRITE  16,106) 

SN  0041  24  CONTINUE 

•:SN  0042  IFIIPRINT.EO.il  RETURN 

C - PRINTOUT  OF  ARRAYS 
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SN  0044 
SN  0045 
SN  0046 
SN  0048 
SN  0049 
SN  0050 
SN  0051 
SN  0052 
SN  0053 
SN  0054 
SN  0055 
SN  0056 
SN  0057 
SN  0058 


>40— ’I'  MlffUAl  MMARGM  CO—UnON 

>1011.7  ~ 

MR  ITI(  *,1001  IA»(1>,ITBAR(I),I-1,MP1) 

100  FORMAT!  * - PROFILES - */lH  ,A8,11FU  .3/|9X,UFU  . 

IFiRHOBAR.EO.l.l  RETURN 

NRITFI6.101I  LASI2) ,|RHOA!I) ,1-1 ,NP1) 

101  F01MAT11H  ,A8,11F11.5/!9X,11F11.5II 
00  30  1-2, N 

J-NP  !♦ 1- 1 

N6( JI-.5*|MG<  JI4N6I J-l I) 

30  CONTINUE 
NGI  11-HGl 

WRITE! 6, 108)  LAB  I  3) , ! NGI I ) , I *1 ,NP1 ) 

1  OR  FORMAT! 1H  ,AB,1P11E11.3/!9X,11E11.3)I 

return 

END 
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